python4oceanographers

Turning ripples into waves

Using Matlab Toolboxes inside IPython

Tidal analysis in python can be quite difficult. During my MatlabTM times I used the t_tide a lot and I end up comparing other packages to t_tide.

The only python alternative one that comes close to t_tide is tappy.

Tappy seems to be a good one, but it was build to be a command line interface, therefore it lacks an interactive mode. Also, it is missing several of t_tide's features, like analyzing a complex vector ($u+iv$).

Because of these shortcomings I still rely on t_tide for my analyses. This post is to show how you can call t_tide toolbox from python using octave.

We will use elevation data from the GLOSS Station 194, loading it in a pandas DataFrame to make our job easier.

In [2]:
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, hour))
    return datetime(year, day, month, hour)

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

obs = read_table('./data/can1998.dtf', names=names, skipinitialspace=True, delim_whitespace=True,
                 index_col='datetime', usecols=range(1, 7), na_fvalues='9.990',
                 parse_dates={'datetime': ['year', 'month', 'day','hour']}, date_parser=date_parser)
obs.head(6)
Out[2]:
elev flag
datetime
1998-01-01 00:00:00 1.20 0
1998-01-01 01:00:00 1.43 0
1998-01-01 02:00:00 1.73 0
1998-01-01 03:00:00 2.03 0
1998-01-01 04:00:00 2.38 0
1998-01-01 05:00:00 2.54 0

6 rows × 2 columns

The metadata explains that there are two flag values, one for bad data and another one for corrected data. A quick plot can help us visualize the flags.

In [3]:
import numpy as np
import matplotlib.pyplot as plt

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

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

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

ax = obs['anomaly'].plot(figsize=(8, 3), alpha=0.5)
ax.set_xlabel('')
_ = ax.set_ylabel(u'Height [m]')
In [4]:
obs['elev'].describe()
Out[4]:
count    8750.000000
mean        1.769914
std         0.393997
min         0.630000
25%         1.480000
50%         1.770000
75%         2.040000
max         3.130000
Name: elev, dtype: float64

The data seems to be OK. Now we can go to the next step.

To call octave we will use IPython's octave extension. The extension relies on the oct2py module to create a two-way connection with octave, so you will need to install octpy.

t_tide needs pwelch, so we need to install octave's signal toolbox and its dependencies as well.

For that, just type the lines below inside an octave session.

pkg install -forge specfun
pkg install -forge control
pkg install -forge general
pkg install -forge signal

We are almost there, the last two steps are a just a few changes in the t_tide toolbox to make it more "octave friendly."

  • In t_tide.m we need to change pwelch's overlap argument from ceil(nx/2) to something between 0, 0.95. I ended up leaving it blank [] (the default).

  • The last case block in t_getconsts.m conflicts with octave internal ctime. Since it is not usually used I just commented it out.

The line below is just to pass the elevation data as a numpy array to oct2py.

In [5]:
elev = obs['anomaly'].values
In [6]:
%load_ext oct2py.ipython
In [7]:
%%octave -i elev -o tidestruc -o pout

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

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

   number of standard constituents used: 59
   Points used: 8759 of 8760
   percent of var residual after lsqfit/var original: 21.11 %
   Greenwich phase computed with nodal corrections applied to amplitude
 and phase relative to center time
   Using nonlinear bootstrapped error estimates
   Generating prediction with nodal corrections, SNR is 2.000000
warning: isstr is obsolete and will be removed from a future version of Octave, please use ischar instead
   percent of var residual after synthesis/var original: 22.25 %
-----------------------------------
date: 01-May-2014
nobs = 8760,  ngood = 8759,  record length (days) = 365.00
start time: 01-Jan-1998
rayleigh criterion = 1.0
Greenwich phase computed with nodal corrections applied to amplitude \n and phase relative to center time

x0= 0.000362, x trend= 0

var(x)= 0.15508   var(xp)= 0.12053   var(xres)= 0.03451
percent var predicted/var original= 77.7 %

     tidal amplitude and phase with 95% CI estimates

tide   freq       amp     amp_err    pha    pha_err     snr
 SSA  0.0002282    0.0198    0.034   321.43   121.26     0.33
 MSM  0.0013098    0.0135    0.031   292.99   140.53     0.19
 MM   0.0015122    0.0209    0.035    43.82   116.16     0.35
 MSF  0.0028219    0.0366    0.039    30.04    66.25     0.86
 MF   0.0030501    0.0330    0.040   266.50    80.54     0.69
 ALP1 0.0343966    0.0012    0.004    97.95   183.29    0.077
 2Q1  0.0357064    0.0010    0.004    63.55   209.97    0.055
 SIG1 0.0359087    0.0002    0.004   141.96   276.36   0.0037
*Q1   0.0372185    0.0267    0.006    65.19    12.94       21
 RHO1 0.0374209    0.0072    0.006    60.58    56.20      1.5
*O1   0.0387307    0.1094    0.006    85.45     2.93    3e+02
 TAU1 0.0389588    0.0007    0.003     4.04   184.25    0.083
 BET1 0.0400404    0.0030    0.005    98.29   132.33     0.35
*NO1  0.0402686    0.0067    0.005    29.19    39.03      2.2
 CHI1 0.0404710    0.0049    0.006    82.72    66.89     0.63
*P1   0.0415526    0.0261    0.005   153.06    10.24       33
*K1   0.0417807    0.0651    0.006   147.85     4.80  1.3e+02
 PHI1 0.0420089    0.0010    0.003    33.20   194.98    0.088
 THE1 0.0430905    0.0082    0.006    32.59    43.08      1.8
 J1   0.0432929    0.0030    0.005   254.65    90.62     0.36
 SO1  0.0446027    0.0047    0.005   225.78    70.47     0.76
 OO1  0.0448308    0.0032    0.006   345.09   125.01     0.26
 UPS1 0.0463430    0.0030    0.007   275.24   136.27     0.19
 OQ2  0.0759749    0.0014    0.004   177.29   180.58     0.16
 EPS2 0.0761773    0.0034    0.004   102.80    93.01     0.69
*2N2  0.0774871    0.0159    0.005   145.60    18.24      9.4
*MU2  0.0776895    0.0258    0.006   156.34    10.71       21
*N2   0.0789992    0.0638    0.005   161.99     4.56  1.5e+02
 NU2  0.0792016    0.0067    0.006   170.25    45.43      1.4
*M2   0.0805114    0.3713    0.005    97.71     0.87  5.2e+03
*MKS2 0.0807396    0.0094    0.006   173.60    43.58      2.1
 LDA2 0.0818212    0.0024    0.004    18.70   116.80     0.35
*L2   0.0820236    0.0322    0.006    60.54    10.92       33
*S2   0.0833333    0.2434    0.006   103.25     1.22  1.9e+03
*K2   0.0835615    0.0773    0.007    91.74     4.37  1.4e+02
 MSN2 0.0848455    0.0026    0.004   267.43   102.01     0.45
 ETA2 0.0850736    0.0051    0.007    56.55    89.91      0.6
*MO3  0.1192421    0.0254    0.008     5.19    15.55      9.4
*M3   0.1207671    0.0678    0.006   232.35     5.06  1.3e+02
*SO3  0.1220640    0.0273    0.008   158.66    16.12       12
*MK3  0.1222921    0.0265    0.006   129.52    15.90       18
*SK3  0.1251141    0.0225    0.007   303.47    19.00      9.3
*MN4  0.1595106    0.0312    0.004   161.91     6.78       78
*M4   0.1610228    0.0567    0.003   213.08     3.06  2.9e+02
 SN4  0.1623326    0.0009    0.003   315.93   166.42     0.11
*MS4  0.1638447    0.0275    0.004   333.76     7.31       61
*MK4  0.1640729    0.0096    0.004   327.38    29.36      4.7
 S4   0.1666667    0.0032    0.003   102.96    70.21      1.1
 SK4  0.1668948    0.0019    0.003     0.91   123.30     0.31
*2MK5 0.2028035    0.0032    0.002   169.26    39.88      2.5
 2SK5 0.2084474    0.0006    0.002   189.00   174.28     0.11
 2MN6 0.2400221    0.0014    0.001   308.47    45.40      1.7
*M6   0.2415342    0.0030    0.001   233.51    21.94        7
*2MS6 0.2443561    0.0061    0.001   254.16    10.65       27
 2MK6 0.2445843    0.0016    0.001   264.05    58.54      1.1
*2SM6 0.2471781    0.0029    0.001   295.81    26.40      6.1
*MSK6 0.2474062    0.0023    0.001   284.64    37.06      2.8
 3MK7 0.2833149    0.0014    0.001    10.38    52.10      1.2
*M8   0.3220456    0.0017    0.001   229.43    34.84      2.9
In [8]:
fig, (ax0, ax1, ax2) = plt.subplots(nrows=3, sharey=True, sharex=True, figsize=(10, 4))

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')

No MatlabTM code was written to produce this post. Our Open Source soul can be at ease ;)

More details regarding t_tide can be found here:

R. Pawlowicz, B. Beardsley, and S. Lentz, "Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE", Computers and Geosciences 28 (2002), 929-937.

In [9]:
HTML(html)
Out[9]:

This post was written as an IPython notebook. It is available for download or as a static html.

Creative Commons License
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

Posted by Filipe Fernandes

Comments