In [1]:
%matplotlib inline
import pandas as pd
import unicodedata
import numpy as np
import os
import glob
import re
import xmltodict
import matplotlib.pyplot as plt
import matplotlib.dates as dates
import matplotlib.ticker as tick
from matplotlib.backends.backend_pdf import PdfPages
import snake
from datetime import datetime
from pylab import rcParams
rcParams['figure.figsize'] = 15, 10
In [2]:
Drive = 'U'
Use g[wellinfo[wellinfo['Well']==wellname]['closest_baro']] instead if you want to match the closest barometer to the data
In [3]:
folder = Drive + ':\\GWP\\Groundwater\\Courthouse_Wash_Transducer_Data'
wellinfofile = Drive + ':\\GWP\\Groundwater\\Courthouse_Wash_Transducer_Data\\wellinfo.csv'
In [4]:
wellinfo = pd.read_csv(wellinfofile,header=0,index_col=0)
wellinfo["G_Elev_m"] = wellinfo["GroundElevation"]/3.2808
wellinfo["Well"] = wellinfo['Well'].apply(lambda x: str(x).lower().strip())
wellinfo['WellID'] = wellinfo.index.values
In [5]:
wellinfo
Out[5]:
In [6]:
manualwls = Drive + ':\\GWP\\Groundwater\\Courthouse_Wash_Transducer_Data\\manual_measurements.csv'
manual = pd.read_csv(manualwls, skiprows=0, parse_dates=1, index_col="DateTime", engine="python")
In [7]:
barofile = folder + '\\*baro*'
barodata = snake.compilation(barofile)
In [8]:
barodata['Level']
barodata.tail(15)
Out[8]:
In [9]:
barodata = snake.hourly_resample(barodata)
In [10]:
infile = Drive + ":\\GWP\\Groundwater\\Courthouse_Wash_Transducer_Data\\"
pathlist = os.path.splitext(infile)[0].split('\\')
barodata.to_csv(pathlist[0] + '\\' + pathlist[1] + '\\' + pathlist[2] + '\\' + pathlist[3] + '\\' + pathlist[4] + '\\' + 'baro' + '.csv')
In [11]:
barodata[['Level']].plot()
Out[11]:
In [12]:
folder = Drive + ":\\GWP\\Groundwater\\Courthouse_Wash_Transducer_Data"
wellinfo = snake.make_files_table(folder, wellinfo)
wellinfo.tail(n=3)
Out[12]:
In [13]:
manual.reset_index(inplace=True)
manual.set_index('WellID',inplace=True)
manual["MeasuredLevel"] = np.nan
manual["Temp"] = np.nan
manual["BaroEfficiencyCorrected"] = np.nan
manual["DeltaLevel"] = np.nan
manual["DriftCorrection"] = np.nan
manual["MeasuredBy"] = np.nan
manual["Tape"] = 1
manual["DTWBelowGroundSurface"] = np.nan
manual["WaterElevation"] = np.nan
#manualrecent["DTWBelowGroundSurface"] = np.nan
manlist= [int(i) for i in manual.index.tolist()]
print manlist
for i in manlist:
try:
manual.loc[i,"DTWBelowCasing"] = manual.loc[i,"MeasuredDTW"]
manual.loc[i,"DTWBelowGroundSurface"] = manual.loc[i,"MeasuredDTW"] - wellinfo.loc[i,"Offset"]
manual.loc[i,"WaterElevation"] = wellinfo.loc[i,'GroundElevation'] - manual.loc[i,"DTWBelowGroundSurface"]
except(KeyError):
pass
outpath = pathlist[0] + '\\' + pathlist[1] + '\\' + pathlist[2] + '\\' + pathlist[3] + '\\' + pathlist[4] + '\\' + 'Manual' + '.csv'
manual.to_csv(outpath, index=True, columns= ["DateTime","MeasuredLevel","Temp","BaroEfficiencyCorrected","DeltaLevel",
"MeasuredDTW","DriftCorrection","DTWBelowCasing","DTWBelowGroundSurface",
"WaterElevation","Tape","MeasuredBy"])
In [14]:
pdf_pages = PdfPages(folder+'wells.pdf')
for i in wellinfo.loc[:,'full_file_name']:
print i
g = snake.imp_new_well(folder+'\\'+i, wellinfo, manual, barodata)
glist = g.columns.tolist()
for j in range(len(glist)):
if 'pw' in glist[j]:
h = glist[j]
y1 = g['WaterElevation'].values
y2 = g[h].values
x1 = g['DateTime'].values
wellname, wellid = snake.getwellid(folder+'\\'+i,wellinfo1)
ylast = wellinfo[wellinfo['WellID']==wellid]['GroundElevation'].values[0] + wellinfo[wellinfo['WellID']==wellid]['Offset'].values[0] - fcl(manual[manual['WellID']== wellid],max(g.index.to_datetime()))[1]
yfirst = wellinfo[wellinfo['WellID']==wellid]['GroundElevation'].values[0] + wellinfo[wellinfo['WellID']==wellid]['Offset'].values[0] - fcl(manual[manual['WellID']== wellid],min(g.index.to_datetime()))[1]
xlast = (fcl(manual[manual['WellID']== wellid],max(g.index.to_datetime()))).name.to_datetime()
xfirst = (fcl(manual[manual['WellID']== wellid],min(g.index.to_datetime()))).name.to_datetime()
x4 = [xfirst,xlast]
y4 = [yfirst,ylast]
fig, ax1 = plt.subplots()
ax1.scatter(x4,y4,color='purple')
ax1.plot(x1,y1,color='blue',label='Water Level Elevation')
ax1.set_ylabel('Water Level Elevation',color='blue')
y_formatter = tick.ScalarFormatter(useOffset=False)
ax1.yaxis.set_major_formatter(y_formatter)
ax2 = ax1.twinx()
ax2.set_ylabel('Barometric Pressure (ft)', color='red')
ax2.plot(x1,y2,color='red',label='Barometric pressure (ft)')
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, loc=3)
plt.title('Well: ' + wellname.title() + ' ' + 'Total Drift = ' + str(g['DriftCorrection'][-1]))
pdf_pages.savefig(fig)
fig.close()
pdf_pages.close()
In [22]:
infile = Drive + ':\\GWP\\Groundwater\\Courthouse_Wash_Transducer_Data\\chwnav 2015_06_12.xle'
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[str(obj['Body_xle']['Ch2_data_header']['Identification']).title()] = f['ch2']
tempunit = (obj['Body_xle']['Ch2_data_header']['Unit'])
if tempunit == 'Deg C' or tempunit == u'\N{DEGREE SIGN}' + u'C':
f[str(obj['Body_xle']['Ch2_data_header']['Identification']).title()] = f['ch2'].convert_objects(convert_numeric=True)
elif tempunit == 'Deg F' or tempunit == u'\N{DEGREE SIGN}' + u'F':
f[str(obj['Body_xle']['Ch2_data_header']['Identification']).title()] = f['ch2'].convert_objects(convert_numeric=True)*0.33456
unit = str(obj['Body_xle']['Ch1_data_header']['Unit']).lower()
if unit == "feet" or unit == "ft":
f[str(obj['Body_xle']['Ch1_data_header']['Identification']).title()] = f['ch1'].convert_objects(convert_numeric=True)
elif unit == "kpa":
f[str(obj['Body_xle']['Ch1_data_header']['Identification']).title()] = f['ch1'].convert_objects(convert_numeric=True)*0.33456
elif unit == "mbar":
f[str(obj['Body_xle']['Ch1_data_header']['Identification']).title()] = f['ch1'].convert_objects(convert_numeric=True)*0.0334552565551
elif unit == "psi":
f[str(obj['Body_xle']['Ch1_data_header']['Identification']).title()] = f['ch1'].convert_objects(convert_numeric=True)*2.306726
elif unit == "m" or unit == "meters":
f[str(obj['Body_xle']['Ch1_data_header']['Identification']).title()] = f['ch1'].convert_objects(convert_numeric=True)*3.28084
else:
f[str(obj['Body_xle']['Ch1_data_header']['Identification']).title()] = f['ch1'].convert_objects(convert_numeric=True)
print "Unknown Units"
# add extension-free file name to dataframe
f['name'] = snake.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[str(obj['Body_xle']['Ch1_data_header']['Identification']).title()] = f[str(obj['Body_xle']['Ch1_data_header']['Identification']).title()].convert_objects(convert_numeric=True)
f[str(obj['Body_xle']['Ch2_data_header']['Identification']).title()] = f[str(obj['Body_xle']['Ch2_data_header']['Identification']).title()].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)
In [23]:
f
Out[23]:
In [42]:
import snake
inputfile = r"U:\GWP\Groundwater\Courthouse_Wash_Transducer_Data\chwnav 2015_06_12.xle"
g = snake.imp_new_well(inputfile, wellinfo, manual, barodata)
glist = g.columns.tolist()
for j in range(len(glist)):
if 'pw' in glist[j]:
h = glist[j]
y1 = g['WaterElevation'].values
y2 = g[h].values
x1 = g['DateTime'].values
wellname, wellid = getwellid(inputfile)
ylast = wellinfo[wellinfo['WellID']==wellid]['GroundElevation'].values[0] + wellinfo[wellinfo['WellID']==wellid]['Offset'].values[0] - fcl(manual[manual['WellID']== wellid],max(g.index.to_datetime()))[1]
yfirst = wellinfo[wellinfo['WellID']==wellid]['GroundElevation'].values[0] + wellinfo[wellinfo['WellID']==wellid]['Offset'].values[0] - fcl(manual[manual['WellID']== wellid],min(g.index.to_datetime()))[1]
xlast = (fcl(manual[manual['WellID']== wellid],max(g.index.to_datetime()))).name.to_datetime()
xfirst = (fcl(manual[manual['WellID']== wellid],min(g.index.to_datetime()))).name.to_datetime()
x4 = [xfirst,xlast]
y4 = [yfirst,ylast]
fig, ax1 = plt.subplots()
ax1.scatter(x4,y4,color='purple')
ax1.plot(x1,y1,color='red')
y_formatter = tick.ScalarFormatter(useOffset=False)
ax1.yaxis.set_major_formatter(y_formatter)
ax2 = ax1.twinx()
ax2.plot(x1,y2,color='blue')
plt.title(getfilename(inputfile)+' '+str(g['DriftCorrection'][-1]))
plt.show()
In [50]:
wellname, wellid = snake.getwellid(infile,wellinfo) #see custom getwellid function
f = snake.new_xle_imp(infile)
bse = int(f.index.to_datetime().minute[0])
try:
bp = wellinfo[wellinfo['Well']==wellname]['BE barologger']
b = snake.hourly_resample(barodata[bp], bse)
b = b.to_frame()
except (KeyError,NameError):
bp = u'Level'
b = snake.hourly_resample(barodata['Level'], bse)
b = b.to_frame()
b['bp'] = b['Level']
b.drop(['Level'],inplace=True,axis=1)
f = snake.hourly_resample(f,bse)
print f
g = pd.merge(f,b,left_index=True,right_index=True,how='inner',sort=True)
print g
g['MeasuredLevel'] = g['Level']
# Remove first and/or last measurements if the transducer was out of the water
## Examine First Value
firstupper = np.mean(g['MeasuredLevel'].diff()[2:31]) + np.std(g['MeasuredLevel'].diff()[2:31])*2.2 # 2.2 std dev.
firstlower = np.mean(g['MeasuredLevel'].diff()[2:31]) - np.std(g['MeasuredLevel'].diff()[2:31])*2.2 # 2.2 std dev.
firstlev = g['MeasuredLevel'].diff()[1:2].values[0] # difference of first two values
## Examine Last Value
lastupper = np.mean(g['MeasuredLevel'].diff()[-31:-2]) + np.std(g['MeasuredLevel'].diff()[-31:-2])*2.2 # 2.2 std dev.
lastlower = np.mean(g['MeasuredLevel'].diff()[-31:-2]) - np.std(g['MeasuredLevel'].diff()[-31:-2])*2.2 # 2.2 std dev.
lastlev = g['MeasuredLevel'].diff()[-2:-1].values[0] # difference of last two values
## drop first value if 2.2 std dev beyond first 30 values
if np.abs(firstlev) > 0.1:
if firstlev > firstupper or firstlev < firstlower:
g.drop(g.index[0],inplace=True)
## drop last value if 2.2 std dev beyond last 30 values
if np.abs(lastlev) > 0.1:
if lastlev > lastupper or lastlev < lastlower:
g.drop(g.index[-1],inplace=True)
glist = f.columns.tolist()
if 'Temperature' in glist:
g['Temp'] = g['Temperature']
g.drop(['Temperature'],inplace=True,axis=1)
elif 'Temp' in glist:
pass
# Get Baro Efficiency
be = wellinfo[wellinfo['WellID']==wellid]['BaroEfficiency']
be = be.iloc[0]
# Barometric Efficiency Correction
g['BaroEfficiencyCorrected'] = g['MeasuredLevel'] - g[bp] + be*g[bp]
In [ ]:
q = {}
for i in wellinfo.loc[:,'full_file_name']:
wellname, wellid = getwellid(folder+'\\'+i)
q[wellname] = imp_new_well(folder+'\\'+i, wellinfo, manual)
In [ ]:
q.names.to_list()
In [ ]:
alldf = ag13a.append(q)
outpath = pathlist[0] + '\\' + pathlist[1] + '\\' + pathlist[2] + '\\' + pathlist[3] + '\\' + pathlist[4] + '\\' + 'all' + '.csv'
alldf.to_csv(outpath, index=False, columns= ["WellID","DateTime","MeasuredLevel","Temp","BaroEfficiencyCorrected","DeltaLevel",
"MeasuredDTW","DriftCorrection","DTWBelowCasing","DTWBelowGroundSurface",
"WaterElevation","Tape","MeasuredBy"])
In [ ]:
dailyavgfiles = r"U:\GWP\Snake Valley Water\Transducer Data\Individual Sites\Daily Averages"
In [ ]:
well = r"U:\GWP\Snake Valley Water\Transducer Data\Individual Sites\Daily Averages\ag13b.xlsx"
In [ ]:
hourly = pd.read_excel(well,'input',parse_dates=True,index_col='DateTime')
hwl = hourly["WaterElevation"].resample('D',how='mean')
hwl = hwl.interpolate(how='Time')
hwl.plot()
In [ ]: