In [19]:
%pylab inline
import os
import sys
import math
import datetime
import pandas as pd
import scipy.linalg as splin
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('%m/%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)
$\hspace{10pt}$ accessed from http://tidesandcurrents.noaa.gov/harcon.html?unit=0&timezone=1&id=8723214&name=Virginia+Key&state=FL
In [20]:
# [name, period, amplitude in m, phase in degrees]
noaa_cons = [['M2',12.421,0.298,256.0],
['S2',12.,0.052,281.3],
['N2',12.658,0.066,239.3],
['K1',23.935,0.031,187.9],
['M4',6.210,0.007,139.6],
['O1',25.819,0.028,218.6],
['M6',4.140,0.011,94.6],
['MK3',8.177,0.004,55.1],
['S4',6,0,0],
['MN4',6.269173724,0,0],
['NU2',12.62600509,0.014,235],
['S6',4,0,0],
['MU2',12.8717576,0.006,265.6],
['2N2',12.90537297,0.009,226.8],
['OO1',22.30608083,0.001,157.1],
['LAM2',12.22177348,0.004,277.2],
['S1',24,0,0],
['M1',24.84120241,0.002,203.2],
['J1',23.09848146,0.002,172.6],
['MM',661.3111655,0,0],
['SSA',4383.076325,0.055,57.8],
['SA',8766.15265,0.083,190.3],
['MSF',354.3670666,0,0],
['MF',327.8599387,0,0],
['RHO',26.72305326,0.001,231.9],
['Q1',26.86835,0.006,231.2],
['T2',12.01644934,0.006,271.6],
['R2',11.98359564,0,282.3],
['2Q1',28.00621204,0.001,249.2],
['P1',24.06588766,0.009,188.3],
['2SM2',11.60695157,0,0],
['M3',8.280400802,0,0],
['L2',12.19162085,0.015,263.3],
['2MK3',8.38630265,0,0],
['K2',11.96723606,0.014,278.4],
['M8',3.105150301,0,0],
['MS4',6.103339275,0.004,177.4]]
nfreq = 8
harmonics = noaa_cons[:nfreq]
In [21]:
def harmonic_analysis(station,pharmonics,t,v):
npfreq = len(pharmonics)
npsamples = v.shape[0]
nf2 = 2 * npfreq
nf21, nf22 = nf2+1, nf2+2
aa = zeros((nf22, nf22), float)
ci = zeros(npfreq, float)
si = zeros(npfreq, float)
x = zeros(nf21, float)
trig = zeros(nf22, float)
amplitude = zeros(npfreq, float)
theta = zeros(npfreq, float)
trig[0] = 1.
for i in xrange(npsamples):
for j in xrange(1,nf2,2):
jh = (j-1) / 2
period = pharmonics[jh][1] * 60. * 60.
if period==0.0: print period
fr = 1. / period
omega = 2. * pi * fr
arc = omega * t[i]
j1 = j + 1
trig[j] = cos(arc)
trig[j1] = sin(arc)
trig[nf22-1] = v[i]
# fill matrix
for i1 in xrange(nf21):
for j1 in xrange(i1,nf22):
aa[i1,j1] = aa[i1,j1] + trig[i1] * trig[j1]
aa[j1,i1] = aa[i1,j1]
# scale diagonal and off diagonals
for i in xrange(nf21):
diagonal = aa[i,i]
aa[i,i] = 1.0
for j in xrange(nf22):
if (j != i):
aa[i,j] = aa[i,j] / diagonal
# solve equations
determinant = splin.det(aa[:nf21,:nf21])
if abs(determinant) < 1.e-20:
print 'singular or ill-conditioned matrix'
print 'determinant = (0}'.format(determinant)
ainv = splin.inv(aa[:nf21,:nf21])
# calculate x
for i in xrange(nf21):
x[i] = 0.
for j in xrange(nf21):
x[i] += ainv[i,j] * aa[j,nf22-1]
# calculate amplitude and phase
for i in xrange(npfreq):
i2 = (i * 2) + 1
ci[i] = x[i2]
si[i] = x[i2+1]
amplitude[i] = sqrt(ci[i] * ci[i] + si[i] * si[i])
ttheta = arctan2(si[i], ci[i]) #* 180. / pi #degrees
if ttheta < 0.:
ttheta += 2. * pi
theta[i] = rad2deg(ttheta)
#--write the data
print station
for idx,[c,period,pa,ptheta] in enumerate(pharmonics):
print ' {0:5s} {1: 7.3f} {2: 7.1f} deg.'.format(c, amplitude[idx], theta[idx])+\
' -- actual {0: 7.1f} deg. difference {1: 7.1f} deg.'.format(ptheta, (ptheta-theta[idx]))
return amplitude, theta
In [22]:
#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='#')
dout = []
for [d,v] in raw:
dt = datetime.datetime.strptime(d, '%m/%d/%Y %H:%M:%S')
dout.append([dt,v])
dout = array(dout)
dout[:,1] -= average(dout[:,1])
return dout
In [23]:
data_dir = os.path.join('data')
tide_file = 'VAKey_6min.csv'
tide = ['tide',0.0,parse_data(os.path.join(data_dir,tide_file))]
time_min = tide[2][0,0]
time_max = tide[2][-1,0]
dt_sec = (tide[2][1,0] - tide[2][0,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)
#--convert from ft to meters
tide[2][:,1] /= 3.28081
#--
de_tide = parse_data(os.path.join(data_dir,'DE_Tide.csv'))
#--plot raw normalized data
t = tide[2]
ax = subplot(1,1,1)
ax.plot(t[:,0], t[:,1], linewidth=2)
ax.plot(de_tide[:,0], de_tide[:,1], linewidth=0.75, color='red')
set_monthsFmt(ax)
ylabel('Elevation, m')
show()
In [24]:
def pair_data(v,base):
a = pd.DataFrame(v[:,1],index=v[:,0])
b = pd.DataFrame(base[:,1],index=base[:,0])
c = pd.concat([a,b], axis=1, join='inner')
t = c.index.to_pydatetime()
pairv = array(zip(t, c.values[:,0]))
pbase = array(zip(t, c.values[:,1]))
return pairv, pbase
def calculate_t(dt):
#--calculate time
total_sec = (dt[-1] - dt[0]).total_seconds()
nsamples = dt.shape[0]
t = linspace(0.,float(total_sec),nsamples)
return t
In [25]:
def create_signal(pharmonics,t):
nsamples = t.shape
v = zeros(nsamples, float)
for idx,[c,period,a,theta,te,dtheta] in enumerate(pharmonics):
period *= 60. * 60.
fr = 1. / period
omega = 2. * pi * fr
theta *= pi / 180.
phase = dtheta * pi / 180.
v += te * a * cos(omega*t - theta - phase)
return v
def signal_error(pharmonics,obsv):
terr = calculate_t(obsv[:,0])
simv = create_signal(pharmonics,terr)
resid = subtract(simv,obsv[:,1])
me = average(resid)
rmse = std(resid)
return me, rmse
In [26]:
#--process tide
station = 'tide'
tide_data = tide[2]
t = calculate_t(tide_data[:,0])
tide_alpha, tide_theta = harmonic_analysis(station,harmonics,t,tide_data[:,1])
#--plot simulated and observed tide
harmonic_results = []
for idx,[c,period,a,tv] in enumerate(harmonics):
harmonic_results.append([c, period, tide_alpha[idx], tide_theta[idx], 1.0, 0.0])
vr = create_signal(harmonic_results,t)
me, rmse = signal_error(harmonic_results,tide_data)
ctxt = '{0} - ME: {1:7.2g} meters, RMSE: {2:7.2g} meters'.format('Tide', me, rmse)
harmonic_results = []
ax = subplot(1,1,1)
ax.plot(tide_data[:,0], vr, linewidth=3.0, color='black', zorder=9)
ax.plot(tide_data[:,0], tide_data[:,1], linewidth=0.75, color='red', zorder=10)
ax.text(0.0, 1.01, ctxt, ha='left', va='bottom', size=8, transform=ax.transAxes)
set_weeksFmt(ax)
#ax.set_xlim(time_min, datetime.date(1999, 2, 1))
show()
ax = subplot(1,1,1)
ax.plot(tide_data[:,1], (vr-tide_data[:,1]), linewidth=0.0, marker='o', markersize=2, color='black')
show()
In [26]:
In [26]:
In [26]: