We will use elevation data from the GLOSS Station 194, loading it in a pandas Series to make our job easier. Check the format of the data below.


In [ ]:
%%bash
more ../../data/gloss_hourly.fmt

In [ ]:
import os
import zipfile
import numpy as np
import numpy.ma as ma
from zipfile import ZipFile
from datetime import datetime
from pandas import Series, date_range


def basename(fname):
    return os.path.splitext(os.path.basename(fname))[0]


def make_date(data):
    date = str(int(data))[:-1]
    return datetime.strptime(date, '%Y%m%d')


def load_gloss(fname):
    data = np.loadtxt(fname, skiprows=1, usecols=range(2, 15))
    elev = ma.masked_equal(data[:, 1:].ravel(), 9999)
    start = make_date(data[0, 0])
    dates = date_range(start, periods=len(elev), freq='H')
    return Series(elev, index=dates)

In [ ]:
series = Series()
with ZipFile('../../data/h281a.zip', 'r') as ziped:
    for fname in sorted(ziped.namelist()):
        if not fname.endswith('/'):  # Skip directories.
            if False:
                with ziped.open(fname) as f:
                    print(f.readline())
            name = basename(fname)
            series = series.append(load_gloss(ziped.open(fname)))

series /= 1e3
series.describe()

In [ ]:
ax = series[::100].plot(figsize=(12, 2))

In [ ]:
obs = series.ix['1998'] - series.ix['1998'].mean()
ax = obs.plot(figsize=(12, 2))
obs.describe()
elev = obs.values

In [ ]:
%load_ext oct2py.ipython

In [ ]:
%%octave -i elev -o tidestruc -o pout

pkg load all
addpath(genpath('../../m_files/t_tide_v1.3beta'))

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

In [ ]:
import matplotlib.pyplot as plt

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

ax0.plot(obs.index, obs, 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-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]