In [2]:

    

%matplotlib inline
import pandas as pd
import numpy as np
import os
import glob
import urllib2
import xmltodict
import matplotlib.pyplot as plt
import statsmodels.tsa.tsatools as tools
from pandas.stats.api import ols
from pylab import rcParams
rcParams['figure.figsize'] = 15, 10


    

In [3]:

    

def readxle(infile):
    def getfilename(path):
        # this function extracts the file name without file path or extension
        return path.split('\\').pop().split('/').pop().rsplit('.', 1)[0]
    # accomodate for files online
    if infile[0:3]=='htt' or infile[0:3]=='ftp':
        response = urllib2.urlopen(infile)
        data = response.read()
        obj = xmltodict.parse(data,encoding="ISO-8859-1")
    else:
        # open text file
        with open(infile) as fd:
            # parse xml
            obj = xmltodict.parse(fd.read(),encoding="ISO-8859-1")
    # navigate through xml to the data
    wellrawdata = obj['Body_xle']['Data']['Log']
    # convert xml data to pandas dataframe
    f = pd.DataFrame(wellrawdata)
    # get header names and apply to the pandas dataframe
    f[obj['Body_xle']['Ch1_data_header']['Identification']] = f['ch1']
    f[obj['Body_xle']['Ch2_data_header']['Identification']] = f['ch2']

    # add extension-free file name to dataframe
    f['name'] = getfilename(infile)
    # combine Date and Time fields into one field
    f['DateTime'] = pd.to_datetime(f.apply(lambda x: x['Date'] + ' ' + x['Time'], 1))
    f[obj['Body_xle']['Ch1_data_header']['Identification']] = f[obj['Body_xle']['Ch1_data_header']['Identification']].convert_objects(convert_numeric=True)
    f = f.reset_index()
    f = f.set_index('DateTime')
    f = f.drop(['Date','Time','@id','ch1','ch2','index','ms'],axis=1)
    return f


    

In [4]:

    

# import relevant data from github
bpfileurl = "https://raw.githubusercontent.com/inkenbrandt/Earth_Tides/master/pw03%20baro%202015-03-04.xle"
wlfileurl = "https://raw.githubusercontent.com/inkenbrandt/Earth_Tides/master/ag13a%202015-03-03.xle"

bp = readxle(bpfileurl)
bp['bp'] = bp['Level']

wl = readxle(wlfileurl)
# unvented transducer provides absolute pressure in feet of water above transducer
wl['wl_and_bp'] = wl['LEVEL']

# merge bp and wl data
data = pd.merge(bp, wl, left_index=True, right_index=True, how='inner')

# calculate vented water level in feet of water above transducer
data['wl'] = data['wl_and_bp']-data['bp']

# drop old fields
data.drop(['Level','LEVEL','Temperature','TEMPERATURE','name_x','name_y','wl_and_bp'], axis=1, inplace=True)
data.dropna(inplace=True)


    

In [5]:

    

# plot the data
x = data.index.to_datetime()
y1 = data['wl']
y2 = data['bp']

fig, ax1 = plt.subplots()
plt.title('Barometric Pressure and Water Level over Time')
ax2 = ax1.twinx()
ax1.plot(x,y1,c='r',label='Water Level')
ax1.set_ylabel('Well Transducer Pressure (ft)', color='r') 
ax2.set_ylabel('Barometric Pressure (ft)', color='b') 
ax2.plot(x,y2,c='b',label='Barometric Pressure')
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, loc=3)


    

    
Out[5]:





<matplotlib.legend.Legend at 0x8bd3950>






    

In [6]:

    

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)')


    

    
Out[6]:





<matplotlib.text.Text at 0x88cf030>






    

In [7]:

    

# clark's method
def clarks(data,bp,wl):
    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

m,b,r = clarks(data, 'bp','wl')
print m
print b
print r


    

    




0.627174234161
-0.955536232102
0.993171472694






    

In [9]:

    

def baro_eff(df,bp,wl,lag=100):
    df.dropna(inplace=True)
    dwl = df[wl].diff().values[1:-1]
    dbp = df[bp].diff().values[1:-1]
    #dwl = df[wl].values[1:-1]
    #dbp = df[bp].values[1:-1]
    df['j_dates'] = df.index.to_julian_date()
    lag_time = df['j_dates'].diff().cumsum().values[1:-1]
    df.drop('j_dates',axis=1,in_place=True)
    # Calculate BP Response Function

    ## create lag matrix for regression
    bpmat = tools.lagmat(dbp, lag, original='in')
    ## transpose matrix to determine required length
    ## run least squared regression
    sqrd = np.linalg.lstsq(bpmat,dwl)
    wlls = sqrd[0]
    cumls = np.cumsum(wlls)
    negcumls = [-1*cumls[i] for i in range(len(cumls))]
    ymod = np.dot(bpmat,wlls)
    
    ## resid gives the residual of the bp
    resid=[(dwl[i] - ymod[i]) for i in range(len(dwl))]
    lag_trim = lag_time[0:len(cumls)]
    return negcumls, cumls, ymod, resid, lag_time, dwl, dbp


    

In [10]:

    

negcumls, cumls, ymod, resid, lag_time, dwl, dbp = baro_eff(data, 'bp','wl',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')


    

    
Out[10]:





<matplotlib.text.Text at 0x10e76a10>






    

In [11]:

    

plt.figure()
x = data.index.to_datetime()[1:]
y1 = resid
y2 = data['wl'][1:]

fig, ax1 = plt.subplots()
plt.title('Corrected Water Level')
ax2 = ax1.twinx()
ax1.plot(x,y1,c='r',label='Residual Water Level')
ax1.set_ylabel('Residual Water Level (ft)', color='r') 
ax2.set_ylabel('Uncorrected Water Level (ft)', color='b') 
ax2.plot(x,y2,c='b',label='Uncorrected Water Level')
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, loc=3)


    

    
Out[11]:





<matplotlib.legend.Legend at 0x12862390>






    






<matplotlib.figure.Figure at 0x10bbe790>