In [1]:
%pylab inline
import os
import sys
from pytides.tide import Tide
from matplotlib import dates as mdates
#--default matplotlib
rcParams['mathtext.default'] = 'regular'
rcParams['legend.fontsize'] = 8
rcParams['axes.labelsize'] = 8
rcParams['xtick.labelsize'] = 8
rcParams['ytick.labelsize'] = 8
rcParams['figure.subplot.wspace'] = 0.35
rcParams['figure.subplot.hspace'] = 0.35
months = mdates.MonthLocator() # every day
weeks = mdates.WeekdayLocator(byweekday=MO, interval=1) #--every week
days = mdates.DayLocator() # every day
hours = mdates.HourLocator() # every hour
monthsFmt = mdates.DateFormatter('%m')
daysFmt = mdates.DateFormatter('%d')
weeksFmt = mdates.DateFormatter('%m/%d')
def set_daysFmt(ax):
ax.xaxis.set_major_locator(days)
ax.xaxis.set_minor_locator(hours)
ax.xaxis.set_major_formatter(daysFmt)
def set_weeksFmt(ax):
ax.xaxis.set_major_locator(weeks)
ax.xaxis.set_minor_locator(days)
ax.xaxis.set_major_formatter(weeksFmt)
def set_monthsFmt(ax):
ax.xaxis.set_major_locator(months)
ax.xaxis.set_minor_locator(days)
ax.xaxis.set_major_formatter(monthsFmt)
In [2]:
#function to parse data and normalize it
def parse_data(fname):
data_dtype = dtype([('datetime','a20'),('head','f4')])
raw = loadtxt(fname, delimiter=',', dtype=data_dtype, comments='#')
dtime = []
dv = []
for [d,v] in raw:
dt = datetime.datetime.strptime(d, '%m/%d/%Y %H:%M:%S')
dtime.append(dt)
dv.append(v)
dtime = array(dtime)
dv = array(dv)
dv -= average(dv)
return dtime, dv
In [3]:
data_dir = os.path.join('data')
vakey_time, vakey_tide = parse_data(os.path.join(data_dir,'VAKey_6min.csv'))
vakey_tide /= 3.28081 #--convert from ft to m
time_min = vakey_time[0]
time_max = vakey_time[-1]
dt_sec = (vakey_time[1] - vakey_time[0]).total_seconds()
#--create time_all
tstep = datetime.timedelta(seconds=dt_sec)
dt = time_min
time_all = []
while dt <= time_max:
time_all.append(dt)
dt += tstep
time_all = array(time_all)
#--deering estates tidal data - already in meters
de_time, de_tide = parse_data(os.path.join(data_dir,'DE_Tide.csv'))
#--plot raw normalized data
ax = subplot(1,1,1)
ax.plot(vakey_time, vakey_tide, linewidth=2)
ax.plot(de_time, de_tide, linewidth=0.75, color='red')
set_weeksFmt(ax)
ylabel('Elevation, m')
show()
In [4]:
total_hours = float((time_all[-1] - time_all[0]).total_seconds()) / (60. * 60.)
nsamples = time_all.shape[0]
hours = linspace(0.,float(total_hours),nsamples)
times = Tide._times(time_all[0], hours)
In [5]:
def ptide_error(obsv,simv):
resid = subtract(simv,obsv)
me = average(resid)
rmse = std(resid)
return me, rmse
In [6]:
##Fit the tidal data to the harmonic model using Pytides
pvakey_tide = Tide.decompose(vakey_tide, t=vakey_time)
##Predict the tides using the Pytides model.
pvakey_prediction = pvakey_tide.at(times)
In [7]:
#--plot predicted data
me, rmse = ptide_error(vakey_tide,pvakey_prediction)
ctxt = '{0} - ME: {1:7.2g} meters, RMSE: {2:7.2g} meters'.format('Tide', me, rmse)
ax = subplot(1,1,1)
ax.plot(vakey_time, vakey_tide, linewidth=0.75, color='red', zorder=10)
ax.plot(vakey_time, pvakey_prediction, linewidth=3, color='black', zorder=9)
ax.text(0.0, 1.01, ctxt, ha='left', va='bottom', size=8, transform=ax.transAxes)
set_weeksFmt(ax)
ylabel('Elevation, m')
show()
ax = subplot(1,1,1)
ax.plot(vakey_tide, (pvakey_prediction-vakey_tide), linewidth=0.0, marker='o', markersize=2, color='black')
show()
In [12]:
for idx,c in enumerate(pvakey_tide.model['constituent']):
print ' {0:10s} {1: 7.3f} {2: 7.1f} deg.'.format(c.name, pvakey_tide.model['amplitude'][idx], pvakey_tide.model['phase'][idx])
In [ ]: