In [ ]:

    

from datetime import datetime
from pandas import read_table, DataFrame

def date_parser(year, day, month, hour):
    year, day, month, hour = map(int, (year, day, month, float(hour)))
    return datetime(year, day, month, hour)

names = ['secs', 'year', 'month', 'day', 'hour', 'elev', 'flag']

obs = read_table('../../data/uba1998.dat', names=names, skipinitialspace=True, delim_whitespace=True,
                 index_col='datetime', usecols=range(1, 8), sep='[][\'" ]+',
                 parse_dates={'datetime': ['year', 'month', 'day','hour']}, date_parser=date_parser)
obs.head(12)


    

In [ ]:

    

import numpy as np
import matplotlib.pyplot as plt

bad = obs['flag'] == 2
suspected = obs['flag'] == 1

fig, ax = plt.subplots(figsize=(14, 4))
ax.plot(obs.index, obs['elev'], alpha=0.5)
ax.plot(obs.index[bad], obs['elev'][bad], 'ro', label='bad')
ax.plot(obs.index[suspected], obs['elev'][suspected], 'go', label='suspected')
ax.legend(numpoints=1, loc='upper right')

# Remove the mean, and interpolate the bad data points.
obs[bad] = np.NaN
obs[suspected] = np.NaN
obs['anomaly'] = obs['elev'] - obs['elev'].mean()
obs['anomaly'] = obs['anomaly'].interpolate()  # Interpolate gaps.

ax = obs['anomaly'].plot(figsize=(14, 4), alpha=0.5)
ax.set_xlabel('')
_ = ax.set_ylabel(u'Height [m]')


    

In [ ]:

    

obs['elev'].describe()


    

In [ ]:

    

elev = obs['anomaly'].values


    

In [ ]:

    

%load_ext oct2py.ipython


    

In [ ]:

    

%%octave -i elev -o tidestruc -o pout

pkg load all
addpath(genpath('../../m_files'))

[tidestruc, pout] = t_tide(elev, 'interval', 1, 'latitude', -25, 'starttime', [1998, 1, 1, 0]);


    

In [ ]:

    

fig, (ax0, ax1, ax2) = plt.subplots(nrows=3, sharey=True, sharex=True, figsize=(14, 6))

ax0.plot(obs.index, obs['anomaly'], label=u'Observations')
ax0.legend(numpoints=1, loc='lower right')

ax1.plot(obs.index, pout.squeeze(), alpha=0.5, label=u'Prediction')
ax1.legend(numpoints=1, loc='lower right')

ax2.plot(obs.index, obs['anomaly']-pout.squeeze(), alpha=0.5, label=u'Residue')
_ = ax2.legend(numpoints=1, loc='lower right')


    

In [ ]:

    

(tidestruc['tidecon'][16, 0] + tidestruc['tidecon'][10, 0]) / (tidestruc['tidecon'][29, 0] + tidestruc['tidecon'][33, 0])


    

In [ ]:

    

tidestruc['name'][16], tidestruc['name'][10], tidestruc['name'][29], tidestruc['name'][33]