Come prevedere le maree usando le costanti armoniche

Question

Sto cercando di capire come creare una tabella di predizione della marea usando le costanti armoniche.

ho usato le costanti armoniche da Bridgeport (http://tidesandcurrents.noaa.gov/data_menu.shtml?stn=8467150%20Bridgeport,%20CT&type=Harmonic%20Constituents)

e sommati i componenti di marea che utilizzano questo script python -

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 

che stampa un'altezza all'ora per 'oggi'. Le altezze risultanti sono in circa la gamma che mi aspetterei, da -0,5 a 7,5 piedi, ma non sono corrette per la data.

Sono sulla buona strada? Come posso determinare l'epoca della marea? Nell'esempio di Doodsen su wikipedia (http://en.wikipedia.org/wiki/Arthur_Thomas_Doodson) hanno usato 0 per ottenere un risultato per il 1 ° settembre 1991. Ma hanno valori armonici diversi rispetto a quelli pubblicati correnti, e quella data non sembra funzionare per me.

Ecco il contenuto del mio file bridgeport.txt -

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

Se si sta utilizzando C Un approccio potrebbe essere per utilizzare libTCD http://www.flaterco.com/xtide/libtcd.html per memorizzare i dati costituenti, questo fornisce un buon framework per l'accesso ai dati.

Per eseguire i calcoli è necessario regolare l'epoca in modo che sia l'inizio dell'anno corrente. Questo codice utilizza le funzioni di libTCD e si basa sugli algoritmi utilizzati 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; 
} 

quindi calcolare la marea per l'offset dall'inizio dell'anno:

- (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 2

Il National Tidal Datum Epoch non è legato alla Unix epoch; è un riferimento di grandezza del livello medio di marea in quel periodo. I dati di Bridgeport contengono una grandezza e una fase per ciascun componente, quindi forniscono la frequenza del componente nell'ultima colonna. La fase è relativa al GMT (GMT è 0 gradi), quindi specificare il tempo in ore in GMT o compensarlo in base al fuso orario. L'ultima colonna è una funzione del componente, non della posizione, quindi quella colonna è la stessa per qualsiasi posizione. Quindi il problema è sommare sinusoidi di frequenze diverse, come l'analisi di Fourier. Molto probabilmente la funzione risultante non ha un periodo di 24 ore, poiché i componenti hanno tutti periodi diversi (non di 24 ore). Ho usato this come riferimento.

Nel codice seguente, ho incluso alcuni argomenti facoltativi per aggiungere alcuni modificatori di offset (come il tuo 3.16999, non so da dove provenga). Ho separato il codice di lettura e analisi dal codice di calcolo. La funzione restituita da get_tides contiene i dati misurati in un closure. Non devi farlo in quel modo e potresti usare invece una classe.

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

E l'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: per una data specifica, oppure il tempo:

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

Answer 3

Le "Fasi sono in gradi, riferiti a GMT "sul tuo link significa che le maree di Bridgeport sono gradi" di fase "fuori sincrono dalle maree teoriche di" equilibrio "al meridiano di Greenwich.

Avrete anche bisogno delle maree di equilibrio in funzione del tempo. Un metodo consisterebbe nel copiare le costanti dalla colonna 2010 di http://coastalengineeringmanual.tpub.com/Part-II-Chap5/Part-II-Chap50024.htm e utilizzare 2010-01-01T00: 00: 00 come origine temporale/epoca/zero.

Inoltre, esiste un fattore nodale che modula un numero di componenti a seconda della precessione nodale della luna.

Invece di copiare dalla tabella, è possibile modificare ed eseguire tide_fac.f programma da http://adcirc.org/home/related-software/adcirc-utility-programs/ per produrre un insieme di fattori nodali e costituenti di marea di equilibrio al proprio dato orario/epoca/zero.

Ecco l'output di un programma modificato tide_fac.f per produrre gli argomenti di fase nodali e di marea di equilibrio per 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

Ho appena scritto una piccola libreria Python chiamato Pytides per l'analisi e la previsione delle maree. È ancora nelle prime fasi, ma puoi ottenere ciò che vuoi abbastanza facilmente usandolo. Ho scritto an example su come farlo.

Come diceva Dave, il motivo per cui il tuo metodo non funziona è che le fasi NOAA sono pubblicate relativamente agli argomenti di equilibrio dei costituenti, quindi hai bisogno di un modo per elaborare questi argomenti di equilibrio (suggerisco di far fare Pytides per tu!).