Commit 2ad49119 authored by Michael Rudolf's avatar Michael Rudolf
Browse files

automatic eol-change due to linux

parent 640eb1b6
# ignore the cached variables from running the script
Tests/
.ipynb_checkpoints/
Examples/
# ignore pictures and pdfs
*.pdf
*.png
*.jpg
*.eps
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
*$py.class
# C extensions
*.so
# Distribution / packaging
.Python
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
wheels/
share/python-wheels/
*.egg-info/
.installed.cfg
*.egg
MANIFEST
# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec
# Installer logs
pip-log.txt
pip-delete-this-directory.txt
# Unit test / coverage reports
htmlcov/
.tox/
.nox/
.coverage
.coverage.*
.cache
nosetests.xml
coverage.xml
*.cover
.hypothesis/
.pytest_cache/
# Translations
*.mo
*.pot
# Django stuff:
*.log
local_settings.py
db.sqlite3
# Flask stuff:
instance/
.webassets-cache
# Scrapy stuff:
.scrapy
# Sphinx documentation
docs/_build/
# PyBuilder
target/
# Jupyter Notebook
.ipynb_checkpoints
# IPython
profile_default/
ipython_config.py
# pyenv
.python-version
# celery beat schedule file
celerybeat-schedule
# SageMath parsed files
*.sage.py
# Environments
.env
.venv
env/
venv/
ENV/
env.bak/
venv.bak/
# Spyder project settings
.spyderproject
.spyproject
# Rope project settings
.ropeproject
# mkdocs documentation
/site
# mypy
.mypy_cache/
.dmypy.json
dmypy.json
# Pyre type checker
.pyre/
# ignore the cached variables from running the script
Tests/
.ipynb_checkpoints/
Examples/
# ignore pictures and pdfs
*.pdf
*.png
*.jpg
*.eps
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
*$py.class
# C extensions
*.so
# Distribution / packaging
.Python
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
wheels/
share/python-wheels/
*.egg-info/
.installed.cfg
*.egg
MANIFEST
# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec
# Installer logs
pip-log.txt
pip-delete-this-directory.txt
# Unit test / coverage reports
htmlcov/
.tox/
.nox/
.coverage
.coverage.*
.cache
nosetests.xml
coverage.xml
*.cover
.hypothesis/
.pytest_cache/
# Translations
*.mo
*.pot
# Django stuff:
*.log
local_settings.py
db.sqlite3
# Flask stuff:
instance/
.webassets-cache
# Scrapy stuff:
.scrapy
# Sphinx documentation
docs/_build/
# PyBuilder
target/
# Jupyter Notebook
.ipynb_checkpoints
# IPython
profile_default/
ipython_config.py
# pyenv
.python-version
# celery beat schedule file
celerybeat-schedule
# SageMath parsed files
*.sage.py
# Environments
.env
.venv
env/
venv/
ENV/
env.bak/
venv.bak/
# Spyder project settings
.spyderproject
.spyproject
# Rope project settings
.ropeproject
# mkdocs documentation
/site
# mypy
.mypy_cache/
.dmypy.json
dmypy.json
# Pyre type checker
.pyre/
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 13 15:13:56 2018
@author: Warsitzka, Michael
This script can be used to convert tdms files of the velocity stepping test into text files and to plot the vst-data.
"""
#%%=============IMPORT=====================================================
import os
import numpy as np
from nptdms import TdmsFile
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from os.path import splitext
import itertools
import csv
#%%==============PARAMETERS=================================================
A = 0.022619 # area of shear zone (= surface area of lid) (m^2)
li = 0.0776 # inner lever (center cell to center of shear zone)
lo = 0.1250 # outer lever (center cell to hinge point)
le= 1.611655493 # load point level
v= 30 # shear velocity (mm/min)
#g = 9.81 # acceleration due to gravity
normalstress= 2000 # normal stress in Pa, only for .txt or asc. files
#%%==================NAMES=================================================
path_in = 'input/' #define path of input data
path_out = 'output/' #define path of input data
curve_fit= 'logx' #define type of curve fit for shear velocity-shear stress curve:
#>>'lin' = linear:
# [slope, intercept] = logfit(x,y,'linear')
# yApprox = (intercept)+(slope)*x
#>> 'logx' = logarithmic:
# [slope, intercept] = logfit(x,y,'logx')
# yApprox = (intercept)+(slope)*log10(x)
#>>'logy' = exponential:
# [slope, intercept] = logfit(x,y,'logy')
# yApprox = (10^intercept)*(10^slope).^x
#>>'loglog' = powerlaw:
# [slope, intercept] = logfit(x,y,'loglog')
# yApprox = (10^intercept)*x.^(slope)
#%%==========================READ & CONVERT DATA==============================================
files = [i for i in os.listdir(path_in) \
if i.endswith(".tdms") and \
('VST' in i or 'vst' in i)]
n = len(files)
names = [splitext(files[i])[0] for i in range(0,n)]
dfvst_raw = np.empty(n, dtype=object)
dfvst = np.empty(n, dtype=object)
i=0
for i in range(0,n):
tdms = TdmsFile(path_in + files[i])
groups = tdms.groups() # get group names
channels = tdms.group_channels(groups[0]) #get channel names
inc = list((tdms.object(groups[0], 'Velocity').properties).values())[2] # time increment between each point
dfvst_raw[i] = tdms.object(groups[0]).as_dataframe() #load tdms as data frame
shearforce = dfvst_raw[i]['Shear Force'] #velocity in [N]
normalforce = dfvst_raw[i]['Normal Force'] #velocity in [N]
liddispl = dfvst_raw[i]['Lid Displacement'] #velocity in [mm]
vel = dfvst_raw[i]['Velocity'] # shear velocity in [mm/s]
if filter(lambda x: x < 0, vel) == True:
vel[vel < 0] = nan
time = np.arange(0, len(vel)*inc, inc) # time in [s]
displ = (np.cumsum(inc*vel)) #cumulative displacement in [mm]
sigma = normalforce/A #calculate normal stress
tau = (shearforce*lo)/(li*A) #calculate shear stress
fric = (tau/sigma) #friction []
dil = -(liddispl-liddispl[0]) # calculation dilation
dfvst[i]= pd.DataFrame.from_dict({'time': time,
'displacement': displ,
'velocity': vel,
'normalstress': sigma,
'shearstress': tau,
'friction': fric,
'dilation':dil}, orient='index').transpose()
#%%==========FIT SHEAR LOAD CURVE============================================
if curve_fit== 'lin': #'linear'
logvel=vel
logfric=fric
elif curve_fit== 'logx': #'logarithmic'
logvel=np.log10(vel)
logfric=fric
elif curve_fit== 'logy': #'exponential'
logfric=np.log10(fric)
logvel=vel
elif curve_fit== 'loglog': #'powerlaw'
logvel=np.log10(vel)
logfric=np.log10(fric)
# linear fit and evaluation:
slope, intercept,r,p,std = stats.linregress(logvel,logfric)
fitfric= np.polyval((slope,intercept),logvel)
# Standard deviations:
m=len(logvel)
s = np.sum((fric-intercept-slope*logvel)**2)/(m-2)
std_slope = np.sqrt(m*s/(m*np.sum(logvel**2)-np.sum(logvel)**2)) #standard error Coef. friction
std_intercept = np.sqrt(np.sum(logvel**2) * s/ (m*np.sum(logvel**2)-np.sum(logvel)**2)) #standard error cohesion
# rescale the approximation results for plotting
if curve_fit == 'lin':
name_fit='y = '+str(round(slope,3))+'*x +'+str(round(intercept,3))
elif curve_fit=='logx':
name_fit = r'y = ({:.4f}$\pm${:.4f})log10(x) + ({:.4f}$\pm${:.4f})'.format(slope,
std_slope,intercept,std_intercept)
if curve_fit=='logy':
fitfric=10**fitfric
logfric=10**logfric
name_fit = 'y = (10^('+str(round(intercept,3))+')*10^('+str(round(slope,3))+')^x'
elif curve_fit=='loglog':
ran=10**ran
fitfric=10**fitfric
logfric=10**logfric
name_fit= 'y = 10^('+str(round(intercept,3))+')*x^('+str(round(slope,3))+')'
# Calculate MSE and R2
MSE = np.mean((logfric-fitfric)**2) # mean squared error.
Covy = np.cov(fitfric, fric) # = cov(estimated values, experimental values)
R2 = (Covy[1]**2)/(np.var(fitfric)*np.var(fric)) #variance
#%%=======================SAVE VST DATA TO TXT FILE=================================================
# data_vst = pd.DataFrame({'Time [s]': time, 'Shear velocity [mm/s]': vel, 'Normal force [N]': normalforce, 'Shear force [N]': shearforce})
data_vst = np.array((time, vel, normalforce, shearforce))
#data_vst.to_csv(path+names[i] + '.txt',index=None, sep='\t', mode='a')
header_vst = 'Time [s]'+'\t' 'Shear velocity [mm s^-1]'+ '\t'+'Normal force [N]'+'\t'+'Shear force [N]'
with open(path_out + names[i] + '.txt', 'w') as f:
f.write(header_vst + '\n')
for row in data_vst.T:
csv.writer(f, delimiter='\t', lineterminator='\n').writerow(row)
f.close
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 13 15:13:56 2018
@author: Warsitzka, Michael
This script can be used to convert tdms files of the velocity stepping test into text files and to plot the vst-data.
"""
#%%=============IMPORT=====================================================
import os
import numpy as np
from nptdms import TdmsFile
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from os.path import splitext
import itertools
import csv
#%%==============PARAMETERS=================================================
A = 0.022619 # area of shear zone (= surface area of lid) (m^2)
li = 0.0776 # inner lever (center cell to center of shear zone)
lo = 0.1250 # outer lever (center cell to hinge point)
le= 1.611655493 # load point level
v= 30 # shear velocity (mm/min)
#g = 9.81 # acceleration due to gravity
normalstress= 2000 # normal stress in Pa, only for .txt or asc. files
#%%==================NAMES=================================================
path_in = 'input/' #define path of input data
path_out = 'output/' #define path of input data
curve_fit= 'logx' #define type of curve fit for shear velocity-shear stress curve:
#>>'lin' = linear:
# [slope, intercept] = logfit(x,y,'linear')
# yApprox = (intercept)+(slope)*x
#>> 'logx' = logarithmic:
# [slope, intercept] = logfit(x,y,'logx')
# yApprox = (intercept)+(slope)*log10(x)
#>>'logy' = exponential:
# [slope, intercept] = logfit(x,y,'logy')
# yApprox = (10^intercept)*(10^slope).^x
#>>'loglog' = powerlaw:
# [slope, intercept] = logfit(x,y,'loglog')
# yApprox = (10^intercept)*x.^(slope)
#%%==========================READ & CONVERT DATA==============================================
files = [i for i in os.listdir(path_in) \
if i.endswith(".tdms") and \
('VST' in i or 'vst' in i)]
n = len(files)
names = [splitext(files[i])[0] for i in range(0,n)]
dfvst_raw = np.empty(n, dtype=object)
dfvst = np.empty(n, dtype=object)
i=0
for i in range(0,n):
tdms = TdmsFile(path_in + files[i])
groups = tdms.groups() # get group names
channels = tdms.group_channels(groups[0]) #get channel names
inc = list((tdms.object(groups[0], 'Velocity').properties).values())[2] # time increment between each point
dfvst_raw[i] = tdms.object(groups[0]).as_dataframe() #load tdms as data frame
shearforce = dfvst_raw[i]['Shear Force'] #velocity in [N]
normalforce = dfvst_raw[i]['Normal Force'] #velocity in [N]
liddispl = dfvst_raw[i]['Lid Displacement'] #velocity in [mm]
vel = dfvst_raw[i]['Velocity'] # shear velocity in [mm/s]
if filter(lambda x: x < 0, vel) == True:
vel[vel < 0] = nan
time = np.arange(0, len(vel)*inc, inc) # time in [s]
displ = (np.cumsum(inc*vel)) #cumulative displacement in [mm]
sigma = normalforce/A #calculate normal stress
tau = (shearforce*lo)/(li*A) #calculate shear stress
fric = (tau/sigma) #friction []
dil = -(liddispl-liddispl[0]) # calculation dilation
dfvst[i]= pd.DataFrame.from_dict({'time': time,
'displacement': displ,
'velocity': vel,
'normalstress': sigma,
'shearstress': tau,
'friction': fric,
'dilation':dil}, orient='index').transpose()
#%%==========FIT SHEAR LOAD CURVE============================================
if curve_fit== 'lin': #'linear'
logvel=vel
logfric=fric
elif curve_fit== 'logx': #'logarithmic'
logvel=np.log10(vel)
logfric=fric
elif curve_fit== 'logy': #'exponential'
logfric=np.log10(fric)
logvel=vel
elif curve_fit== 'loglog': #'powerlaw'
logvel=np.log10(vel)
logfric=np.log10(fric)
# linear fit and evaluation:
slope, intercept,r,p,std = stats.linregress(logvel,logfric)
fitfric= np.polyval((slope,intercept),logvel)
# Standard deviations:
m=len(logvel)
s = np.sum((fric-intercept-slope*logvel)**2)/(m-2)
std_slope = np.sqrt(m*s/(m*np.sum(logvel**2)-np.sum(logvel)**2)) #standard error Coef. friction
std_intercept = np.sqrt(np.sum(logvel**2) * s/ (m*np.sum(logvel**2)-np.sum(logvel)**2)) #standard error cohesion
# rescale the approximation results for plotting
if curve_fit == 'lin':
name_fit='y = '+str(round(slope,3))+'*x +'+str(round(intercept,3))
elif curve_fit=='logx':
name_fit = r'y = ({:.4f}$\pm${:.4f})log10(x) + ({:.4f}$\pm${:.4f})'.format(slope,
std_slope,intercept,std_intercept)
if curve_fit=='logy':
fitfric=10**fitfric
logfric=10**logfric
name_fit = 'y = (10^('+str(round(intercept,3))+')*10^('+str(round(slope,3))+')^x'
elif curve_fit=='loglog':
ran=10**ran
fitfric=10**fitfric
logfric=10**logfric
name_fit= 'y = 10^('+str(round(intercept,3))+')*x^('+str(round(slope,3))+')'
# Calculate MSE and R2
MSE = np.mean((logfric-fitfric)**2) # mean squared error.
Covy = np.cov(fitfric, fric) # = cov(estimated values, experimental values)
R2 = (Covy[1]**2)/(np.var(fitfric)*np.var(fric)) #variance
#%%=======================SAVE VST DATA TO TXT FILE=================================================
# data_vst = pd.DataFrame({'Time [s]': time, 'Shear velocity [mm/s]': vel, 'Normal force [N]': normalforce, 'Shear force [N]': shearforce})
data_vst = np.array((time, vel, normalforce, shearforce))
#data_vst.to_csv(path+names[i] + '.txt',index=None, sep='\t', mode='a')
header_vst = 'Time [s]'+'\t' 'Shear velocity [mm s^-1]'+ '\t'+'Normal force [N]'+'\t'+'Shear force [N]'
with open(path_out + names[i] + '.txt', 'w') as f:
f.write(header_vst + '\n')
for row in data_vst.T:
csv.writer(f, delimiter='\t', lineterminator='\n').writerow(row)
f.close
# -*- coding: utf-8 -*-
"""
Created on Tue Sep 25 14:54:43 2018
@author: analab
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import csv
import pylab
from scipy import stats
import os
import fnmatch
import glob, shutil
import pickle
from os.path import isfile, join, basename, splitext
from nptdms import TdmsFile
import collections
import itertools
import codecs
#%%==========================PARAMETERS & VARIABLES==============================================
A = 0.022619 # area of shear zone (= surface area of lid) (m^2)
li = 0.0776 # inner lever (center cell to center of shear zone)
lo = 0.1250 # outer lever (center cell to hinge point)
le= 1.611655493 #load point level
v= 30 #shear velocity (mm/min)
g = 9.81 # acceleration due to gravity
prec = 125 # rounding precision for normal load in file names
velnoise = 1e-2
stressnoise = 0.025
smoothwindow = 51 # define window for smooth shear stress curve
projectname = '324-01_QuartzneusandBern_basal'
TSname= projectname+'_ts'
path_in = 'input/'
path_out = 'output/'
#%%==========================READ & CONVERT DATA==============================================
for file in os.listdir(path_in):
files,n = glob.glob1(path_in,"*.asc"),len(glob.glob1(path_in,'*.asc'))
DF = np.empty(n, dtype=object)
len_m = np.zeros(n)
Sigma_avg = np.zeros(m)
for i in range(0,len(files)):
if files[i].endswith('.asc'): # loop over all asc files
name = splitext(files[i])[0] # extract names of asc files
with codecs.open(path_in + files[i], encoding='utf-8-sig') as f:
data_load = np.loadtxt(f) # load file for further calculations
# write data to new varables with meaningfull names
time = data_load[:,0]
shearload = data_load[:,1] #in kg
velocity = data_load[:,2]
dilation = data_load[:,3]
normalload = data_load[:,4] #in kg
# Corrections
if time[0] >= 0:
time = time - time[0] # set starting time to 0
else:
time = time
dilation = -(dilation-dilation[0]) # set initial dilation to 0
shearforce = shearload *g #in N
normalforce = normalload *g #in N
sigma = normalforce/A # kg to normal stress in [Pa]
tau = (shearforce*lo)/(li*A) # kg to shear stress in [Pa]
velocity = velocity/60 # transform to [mm/s]
velocity[velocity < 0] = 0
inc = time[1]-time[0] # timesteps, first one = 0
#incdisp = np.diff(time)*velocity[:-1] # incremental displacement in mm
displacement = np.cumsum(inc*velocity[:-1]) # cummulative displacement in mm
Sigma_avg[i] = np.average(sigma[sigma>0]) # only for the file name
DF[i]= pd.DataFrame.from_dict({'time': time,