How to predict tides using harmonic constants

Question

I'm trying to understand how to create a tide prediction table using the harmonic constants.

I used the harmonic constants from Bridgeport (http://tidesandcurrents.noaa.gov/data_menu.shtml?stn=8467150%20Bridgeport,%20CT&type=Harmonic%20Constituents)

And summed the tidal components using this python script -

import math
import time

tidalepoch = 0
epoch = time.mktime(time.gmtime()) - tidalepoch
f = open('bridgeport.txt', 'r')
M_PI = 3.14159
lines = f.readlines()
t = epoch - 24 * 3600
i = -24
while t < epoch:
  height = 0
  for line in lines:
    x = line.split()
    A = float(x[2]) # amplitude
    B = float(x[3]) # phase
    B *=  M_PI / 180.0
    C = float(x[4]) # speed
    C *= M_PI / 648000

    # h = R cost (wt - phi)
    height += A * math.cos(C * t - B)

  print str(i) + " " + str(height + 3.61999)
  i += 1
  t += 3600

That prints one height per hour for 'today'. The resulting heights are in about the range I would expect, -0.5 to 7.5 feet, but aren't correct for the date.

Am I on the right track? How do I determine the tidal epoch? In the Doodsen example on wikipedia (http://en.wikipedia.org/wiki/Arthur_Thomas_Doodson) they used 0 to get a result for Sept 1, 1991. But they have different harmonic values than the current posted ones, and that date does not seem to work for me.

Here is the content of my bridgeport.txt file -

1 M2                       3.251 109.6  28.9841042 
2 S2                       0.515 135.9  30.0000000 
3 N2                       0.656  87.6  28.4397295 
4 K1                       0.318 191.6  15.0410686 
5 M4                       0.039 127.4  57.9682084 
6 O1                       0.210 219.5  13.9430356 
7 M6                       0.044 353.9  86.9523127 
8 MK3                      0.023 198.8  44.0251729 
9 S4                       0.000   0.0  60.0000000 
10 MN4                      0.024  97.2  57.4238337 
11 NU2                      0.148  89.8  28.5125831 
12 S6                       0.000   0.0  90.0000000 
13 MU2                      0.000   0.0  27.9682084 
14 2N2                      0.077  65.6  27.8953548 
15 OO1                      0.017 228.7  16.1391017 
16 LAM2                     0.068 131.1  29.4556253 
17 S1                       0.031 175.5  15.0000000 
18 M1                       0.024 264.4  14.4966939 
19 J1                       0.021 237.0  15.5854433 
20 MM                       0.000   0.0   0.5443747 
21 SSA                      0.072  61.2   0.0821373 
22 SA                       0.207 132.0   0.0410686 
23 MSF                      0.000   0.0   1.0158958 
24 MF                       0.000   0.0   1.0980331 
25 RHO                      0.015 258.1  13.4715145 
26 Q1                       0.059 205.7  13.3986609 
27 T2                       0.054 106.4  29.9589333 
28 R2                       0.004 136.9  30.0410667 
29 2Q1                      0.014 238.8  12.8542862 
30 P1                       0.098 204.1  14.9589314 
31 2SM2                     0.000   0.0  31.0158958 
32 M3                       0.012 200.1  43.4761563 
33 L2                       0.162 134.1  29.5284789 
34 2MK3                     0.015 203.7  42.9271398 
35 K2                       0.150 134.7  30.0821373 
36 M8                       0.000   0.0 115.9364166 
37 MS4                      0.000   0.0  58.9841042

Answer 1

I've just written a small Python library called Pytides for tidal analysis and prediction. It's still in the early stages, but you can achieve what you want fairly easily using it. I've written up an example of how to do so.

As Dave mentioned, the reason your method isn't working is that NOAA phases are published relative to the equilibrium arguments of the constituents, so you need some way of working out these equilibrium arguments (I suggest letting Pytides do it for you!).

Answer 2

If you're using C one approach would be to use libTCD http://www.flaterco.com/xtide/libtcd.html to store the constituent data, this provides a good framework for accessing the data.

To do the calculations you need to adjust the epoch to be the start of the current year. This code uses functions from libTCD, and is based on algorithms used in xTide.

TIDE_RECORD record;
int readStatus = read_tide_record(self.stationID, &record);
int epochStartSeconds = staringSecondForYear(year);

for (unsigned constituentNumber=0; constituentNumber<constituentCount; constituentNumber++) {
    float constituentAmplitude = record.amplitude[constituentNumber];
    float nodeFactor = get_node_factor(constituentNumber, year);
    amplitudes[constituentNumber] =  constituentAmplitude * nodeFactor;

    float constituentEpoch =  deg2rad(-record.epoch[constituentNumber]);
    float equilibrium = deg2rad(get_equilibrium(constituentNumber, year));
    phases[constituentNumber] = constituentEpoch + equilibrium;
}

Then calculate the tide for the offset since the start of the year:

- (float)getHeightForTimeSince1970:(long)date {
  //calculate the time to use for this calculation
  int secondsSinceEpoch = date - epochStartSeconds;

  float height = 0;
  for(int i = 0; i < constituentCount; i++) {
    float degreesPerHour = get_speed(i);
    float radiansPerSecond = degreesPerHour * M_PI / 648000.0;
    height += amplitudes[i] * cos(radiansPerSecond * secondsSinceEpoch + phases[i]);
  }

  height += record.datum_offset;
  return height;
}

Answer 3

The "Phases are in degrees, referenced to GMT" on your link means that the Bridgeport tides are "phase" degrees out of sync from the theoretical "equilibrium" tides at the Greenwich Meridian.

You will also need the equilibrium tides as a function of time. One method would be to copy the constants from the 2010 column of http://coastalengineeringmanual.tpub.com/Part-II-Chap5/Part-II-Chap50024.htm and use 2010-01-01T00:00:00 as your time datum/epoch/zero.

Additionally, there is a nodal factor that modulates a number of the constituents depending on the nodal precession of the moon.

Rather than copy from the table, you could edit and run the tide_fac.f program from http://adcirc.org/home/related-software/adcirc-utility-programs/ to produce a set of nodal factors and equilibrium tidal constituents at your own chosen time datum/epoch/zero.

Here's output from a modified tide_fac.f program to produce the nodal and equilibrium tide phase arguments for 2013-01-01T00:00:00:

TIDAL FACTORS STARTING:  HR- 0.00,  DAY-  1,  MONTH-  1  YEAR- 2013

FOR A RUN LASTING   365.00 DAYS


CONST   NODE     EQ ARG (ref GM)
NAME   FACTOR    (DEG) 


M2    1.02718     270.21
S2    1.00000       0.00
N2    1.02718      16.12
K1    0.92336      17.70
M4    1.05510     180.42
O1    0.87509     248.94
M6    1.08377      90.63
MK3   0.94846     287.91
S4    1.00000       0.00
MN4   1.05510     286.33
NU2   1.02718      73.02
S6    1.00000       0.00
MU2   1.02718     178.93
2N2   1.02718     122.03
OO1   0.63334     333.60
LAMB  1.02718     287.40
S1    1.00000     180.00
M1    0.00000     159.37
J1    0.89282     275.36
MM    1.09369     254.09
SSA   1.00000     201.62
SA    1.00000     280.81
MSF   1.02718      91.28   
MF    0.74476     312.33
RHO1  0.87509      51.75   
Q1    0.87509     354.85
T2    1.00000       2.35
R2    1.00000     177.65
2Q1   0.87509     100.76
P1    1.00000     349.19
2SM2  1.02718      89.79
M3    1.04114     225.31
L2    0.00000     348.15
2MK3  0.97424     162.72
K2    0.81662     214.58
M8    1.11323       0.84
MS4   1.02718     270.21

Answer 4

The National Tidal Datum Epoch is not related to the Unix epoch; it's a magnitude reference of the average tide level over that period. The Bridgeport data contains a magnitude and phase for each component, and then gives the frequency of the component in the last column. The phase is relative to GMT (GMT is 0 degrees), so specify the time in hours in GMT, or offset it according to timezone. The last column is a function of the component, not the location, so that column is the same for any location. So the problem is summing sinusoids of different frequencies, like Fourier analysis. The resulting function will most likely not have a 24 hour period, as the components all have different (non-24 hour) periods. I used this as a reference.

In the code below, I included some optional arguments to add some offset modifiers (like your 3.16999, not sure where that is from). I separated the reading and parsing code from the calculation code. The returned function from get_tides contains the measured data in a closure. You don't have to do it that way and could use a class instead.

from __future__ import division
import math
import time

def get_tides(fn):
    """Read tide data from file, return a function that
    gives tide level for specified hour."""
    with open(fn,'r') as fh:
        lines = fh.readlines()

    measured = []
    for line in lines:
        index,cons,ampl,phase,speed = line.split()
        measured.append(tuple(float(x) for x in (ampl,phase,speed)))

    def find_level(hour,total_offset=0,time_offset=0):
        total = total_offset
        for ampl,phase,speed in measured:
            # you could break out the constituents here
            # to see the individual contributions
            cosarg = speed * (hour + time_offset) + phase
            total += ampl * math.cos(cosarg * math.pi/180) # do rad conversion here
        return total

    return find_level


bridgeport = get_tides('bridgeport.txt')

for hour in range(-24,25,1):
    print("%3sh %7.4f" % (hour,bridgeport(hour,total_offset=3.61999)))

And the output:

-24h -0.5488
-23h -1.8043
-22h -1.7085
-21h -0.3378
-20h  1.8647
-19h  4.4101
-18h  6.8374
-17h  8.5997
-16h  9.1818
-15h  8.4168
-14h  6.5658
-13h  4.1003
-12h  1.5669
-11h -0.3936
-10h -1.2009
 -9h -0.6705
 -8h  0.9032
 -7h  3.0316
 -6h  5.2485
 -5h  7.0432
 -4h  7.8633
 -3h  7.4000
 -2h  5.8028
 -1h  3.5317
  0h  1.1223
  1h -0.8425
  2h -1.7748
  3h -1.3620
  4h  0.2196
  5h  2.4871
  6h  4.9635
  7h  7.2015
  8h  8.6781
  9h  8.9524
 10h  7.9593
 11h  6.0188
 12h  3.6032
 13h  1.2440
 14h -0.4444
 15h -0.9554
 16h -0.2106
 17h  1.4306
 18h  3.4834
 19h  5.5091
 20h  7.0246
 21h  7.5394
 22h  6.8482
 23h  5.1795
 24h  2.9995

EDIT: For a specific date, or the time now:

import datetime

epoch = datetime.datetime(1991,9,1,0,0,0)
now = datetime.datetime.utcnow()
hourdiff = (now-epoch).days*24

for hour in range(0,25,1):
    print("%3sh %7.4f" % (hour,bridgeport(hour+hourdiff)))