This notebook provides Python scripts to import, compile, modify, graph, and export Solinst transducer data.
In [ ]:
%matplotlib inline
import pandas as pd
import numpy as np
import os
import sys
import platform
import glob
import re
import xmltodict
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.dates as dates
import matplotlib.ticker as tick
from matplotlib.backends.backend_pdf import PdfPages
import statsmodels.tsa.tsatools as tools
from pandas.stats.api import ols
from datetime import datetime
from pylab import rcParams
rcParams['figure.figsize'] = 15, 10
In [ ]:
print("Operating System " + platform.system() + " " + platform.release())
print("Python Version " + str(sys.version))
print("Pandas Version " + str(pd.__version__))
print("Numpy Version " + str(np.__version__))
print("Matplotlib Version " + str(matplotlib.__version__))
In [ ]:
#rootname = '/media/p/Transcend/PROJECTS/UMAR/Phase_II/Data/RAW/'
rootname = 'E:/PROJECTS/UMAR/Data/RAW/'
In [ ]:
def Scat(data,bp,wl):
data['dwl'] = data[wl].diff()
data['dbp'] = data[bp].diff()
regression = ols(y=data['dwl'], x=data['dbp'])
m = regression.beta.x
b = regression.beta.intercept
r = regression.r2
#r = (regression.beta.r)**2
plt.scatter(y=data['dwl'], x=data['dbp'])
y_reg = [data['dbp'][i]*m+b for i in range(len(data['dbp']))]
plt.plot(data['dbp'],y_reg,
label='Regression: Y = {m:.4f}X + {b:.5}\nr^2 = {r:.4f}\n BE = {be:.2f} '.format(m=m,b=b,r=r,be=m))
plt.legend()
plt.xlabel('Sum of Barometric Pressure Changes (ft)')
plt.ylabel('Sum of Water-Level Changes (ft)')
In [ ]:
# clark's method
def clarks(data,bp,wl):
'''
clarks method
Input dataframe (data) with barometric pressure (bp) and water level (wl) data
Returns slope, intercept, and r squared value'''
data['dwl'] = data[wl].diff()
data['dbp'] = data[bp].diff()
data['beta'] = data['dbp']*data['dwl']
data['Sbp'] = np.abs(data['dbp']).cumsum()
data['Swl'] = data[['dwl','beta']].apply(lambda x: -1*np.abs(x[0]) if x[1]>0 else np.abs(x[0]), axis=1).cumsum()
plt.figure()
plt.plot(data['Sbp'],data['Swl'])
regression = ols(y=data['Swl'], x=data['Sbp'])
m = regression.beta.x
b = regression.beta.intercept
r = regression.r2
y_reg = [data.ix[i,'Sbp']*m+b for i in range(len(data['Sbp']))]
plt.plot(data['Sbp'],y_reg,
label='Regression: Y = {m:.4f}X + {b:.5}\nr^2 = {r:.4f}\n BE = {be:.2f} '.format(m=m,b=b,r=r,be=m))
plt.legend()
plt.xlabel('Sum of Barometric Pressure Changes (ft)')
plt.ylabel('Sum of Water-Level Changes (ft)')
data.drop(['dwl','dbp','Sbp','Swl'], axis=1, inplace=True)
return m,b,r
I always set my transducers to future start to make the tranducer start on the hour. I also allow the Levelogger to take an instantaneous measurement out of water, and zero the transducer out to accomodate for elevation.
First, we must import all of the relevant data. To properly import transducer data, we need:
If we want to calculate water-level elevation, we also need:
In [ ]:
barofile = new_xle_imp(rootname + "baro_2015-07-16.xle")
barofile2 = pd.read_csv(rootname + "UCC.csv",parse_dates=True,index_col='Day',skiprows=14, na_values=['M','S'])
wellfile = new_xle_imp(rootname +"arnold_well_2015-07-16.xle")
wellfile2 = new_xle_imp(rootname +"arnold_well_2015-04-01.xle")
manualfile = pd.read_excel(rootname +"Manual_Readings.xlsx","Arn_Well",index_col="datetime")
In [ ]:
barofile2['ft_water_bp']= barofile2['Sea Level Pressure']*0.0335 - (31.17 - 4806/826 + 7.8) # convert hPa to ft water
barofile2 = barofile2.interpolate(method='time') # fill NA spots
Concatonate the well files so that they are one seamless file.
In [ ]:
wellfile = pd.concat([wellfile,wellfile2])
wellfile.sort_index(inplace=True)
In [ ]:
wellfile.columns
You should always graph raw data to see if there are any tares in the data from users moving the tranducer placement. Sometimes, the transducer is out of the water when it takes a measurement. These points should be removed or adjusted.
In [ ]:
#http://stackoverflow.com/questions/7733693/matplotlib-overlay-plots-with-different-scales
x1 = wellfile.index.to_datetime() #converts pandas dataframe index into datetime format for graph
x2 = barofile.index.to_datetime()
x3 = manualfile.index.to_datetime()
y1 = wellfile['Level']
y2 = barofile['Level']
y3 = manualfile['dtw_ft']
data = [(x1,y1),(x2,y2),(x3,y3)]
fig, ax = plt.subplots()
# Twin the x-axis twice to make independent y-axes.
axes = [ax, ax.twinx(), ax.twinx()]
# Make some space on the right side for the extra y-axis.
fig.subplots_adjust(right=0.75)
# Move the last y-axis spine over to the right by 20% of the width of the axes
axes[-1].spines['right'].set_position(('axes', 1.2))
# To make the border of the right-most axis visible, we need to turn the frame
# on. This hides the other plots, however, so we need to turn its fill off.
axes[-1].set_frame_on(True)
axes[-1].patch.set_visible(False)
# And finally we get to plot things...
colors = ['Green', 'Red', 'Blue']
labels = ['Levelogger Pressure (ft)','Barologger Pressure (ft)','Manual Readings (ft to water)' ]
marks = ['','','o']
linetypes = ['solid','solid','none']
for ax, color, datum, label, mark, linety in zip(axes, colors, data, labels, marks, linetypes):
ax.plot(datum[0],datum[1], marker=mark, linestyle=linety, color=color, label=label)
ax.set_ylabel(label, color=color)
ax.tick_params(axis='y', colors=color)
h1, l1 = axes[0].get_legend_handles_labels()
h2, l2 = axes[1].get_legend_handles_labels()
h3, l3 = axes[2].get_legend_handles_labels()
axes[0].legend(h1+h2+h3, l1+l2+l3, loc=4)
plt.show()
In [ ]:
print range(-10,10)
This tranducer has a jump in the middle of the data caused by adjustments during manual recordings, as well as a jump at the beginning due to the transducer being out of water at the time of measurement.
In [ ]:
wellfile = smoother(wellfile, 'Level', 30, 3)
wellfile = smoother(wellfile, 'Conductivity', 30, 3)
In [ ]:
wellfile = jumpfix(wellfile,'Level',0.1)
wellfile = jumpfix(wellfile,'Conductivity',0.005)
wellfile['Level'].plot()
Solinst transducers are nonvented, meaning that they measure absolute pressure. When they are submerged in a well, they are measuring the pressure of the water and the atmosphere. In most cases, we are only interested in the pressure that the water exerts, so we have to subtract the pressure that the atmosphere is exerting.
In [ ]:
wellbaro = baro_drift_correct(wellfile,barofile,manualfile)
In [ ]:
wellbaro.columns
In [ ]:
wellbaro['WaterElevation'].plot()
plt.vlines('11/4/2014 11:16',wellbaro['WaterElevation'].min(),wellbaro['WaterElevation'].max(),color='green')
In [ ]:
Scat(wellbaro,'abs_feet_above_barologger','WaterElevation')
In [ ]:
s, m, r = clarks(wellbaro,'abs_feet_above_barologger','WaterElevation')
In [ ]:
negcumls, cumls, ymod, resid, lag_time, dwl, dbp = baro_eff(wellbaro,'abs_feet_above_barologger','WaterElevation',100)
plt.figure()
lag_trim = lag_time[0:len(negcumls)]
plt.scatter(lag_trim*24,negcumls, label='b.p. alone')
plt.xlabel('lag (hours)')
plt.ylabel('barometric response')
ymin = wellbaro['WaterElevation'].min()
fig, ax = plt.subplots()
plt.plot(wellbaro.index[1:-1], resid)
plt.text(x='11/3/2014 1:00',y=ymin+2,s='Injection Began',rotation=90,color='green',fontsize=12)
y_formatter = tick.ScalarFormatter(useOffset=False)
ax.yaxis.set_major_formatter(y_formatter)
plt.vlines('11/4/2014 11:16',ymin+3,wellbaro['WaterElevation'].max(),color='green')
print len(resid)
print len(wellbaro.index[1:-1])
In [ ]:
wellbaro['corrwl'] = wellbaro['WaterElevation'] - wellbaro['abs_feet_above_barologger']*1
manualfile['wlelev'] = 4800-manualfile['dtw_ft']
x1 = wellbaro.index.to_datetime()[1:-1] #converts pandas dataframe index into datetime format for graph
x2 = barofile.index.to_datetime()
x3 = manualfile.index.to_datetime()
y1 = resid
y2 = barofile['Level']
y3 = manualfile['wlelev']
data = [(x1,y1),(x2,y2),(x3,y3)]
fig, ax = plt.subplots()
# Twin the x-axis twice to make independent y-axes.
axes = [ax, ax.twinx(), ax.twinx()]
# Make some space on the right side for the extra y-axis.
fig.subplots_adjust(right=0.75)
# Move the last y-axis spine over to the right by 20% of the width of the axes
axes[-1].spines['right'].set_position(('axes', 1.2))
# To make the border of the right-most axis visible, we need to turn the frame
# on. This hides the other plots, however, so we need to turn its fill off.
axes[-1].set_frame_on(True)
axes[-1].patch.set_visible(False)
# And finally we get to plot things...
colors = ['Green', 'Red', 'Blue']
labels = ['Levelogger Pressure (ft)','Barologger Pressure (ft)','Manual Readings (ft to water)' ]
marks = ['','','o']
linetypes = ['solid','solid','none']
y_formatter = tick.ScalarFormatter(useOffset=False)
for ax, color, datum, label, mark, linety in zip(axes, colors, data, labels, marks, linetypes):
ax.plot(datum[0],datum[1], marker=mark, linestyle=linety, color=color, label=label)
ax.set_ylabel(label, color=color)
ax.tick_params(axis='y', colors=color)
ax.yaxis.set_major_formatter(y_formatter)
h1, l1 = axes[0].get_legend_handles_labels()
h2, l2 = axes[1].get_legend_handles_labels()
h3, l3 = axes[2].get_legend_handles_labels()
axes[0].legend(h1+h2+h3, l1+l2+l3, loc=4)
axes[2].set_ylim(4485,4493)
plt.show()
It is best to set Solinst transducers (Leveloggers) to start at the same time and to measure at the same frequency as your Barologger. Sometimes, this does not happen. To solve mismatches in sampling interval, we can resample the barometer data to same base (start time) and frequency as the transducer.
Using the hourly_resample function above, we can resample each transducer dataset.