Prédire les marées avec 37 nombres

Dans l’univers de la météo, les marées sont particulières. À peu près tous les aspects des systèmes météorologiques sont la résultante d’intéractions complexes: formation de nuage, vents ou pluie. Pour ces dimensions météo, des modèles mathématiques sont requis pour produire des prévisions fiables, souvent limitées à quelques jours. Ces modèles s’appuient sur des tonnes de données calculant énormément d’information qui est constamment mise-à-jour.

À l’inverse, le niveau des marées peut être prédit à quelques pouces d’écart pour une année entière et ce, avec très peu d’information. C’est ce qui les rend si spéciales dans l’univers de la météo. Les marées sont aussi régulières qu’une horloge. En comparaison aux autres phénomènes météos, elles sont faciles à prédire.

La manière dont elles sont prédites est intimement liée à leurs causes, c’est-à-dire les forces gravitationnelles de la lune, du soleil et des autres planètes, chacune tirant l’eau (et la terre) vers elle. Selon la position de la terre, l’eau est soit plus proche, soit plus éloignée des astres, ce qui crée une différence dans ces forces gravitationnelles.

L’effet net est la création d’ourlets aux côtés opposés de la terre: la marée haute. Des facteurs locaux ont aussi un impact: la profondeur de l’eau, la position du lieu étudié sur terre, la résonance induite par la géographie et d’éventuels contre-courants. Tous peuvent affecter la vitesse, la hauteur et la forme des marées. Un excellent résumé de la théorie des marées est dans le vidéo ci-dessous (en anglais).

Quand la lune et le soleil sont enlignés, les forces gravitationnelles s’additionnent, ce qui crée une grande marée. Quand la lune et le soleil sont orthogonales (perpendiculaires) par rapport à la terre, les forces gravitationnelles agissent dans des directions différentes et les marées sont plus faibles (marée morte).

Cela dit, aussi pédagogique que la vidéo puisse être, elle demeure une explication générale. Elle n’explique pas comment faire une table des marées. Pour des besoins de prévision, il s’avère que 37 nombres, des données sur le niveau d’eau et un ordinateur de bureau sont tout ce dont on a de besoin pour faire des prévisions décentes. La décision cruciale à prendre pour fins de prévision est si on préfère des prévisions précises avec une courte fenêtre de prévision, ou des prévisions sur de longues périodes (e.g. une année), mais un peu moins précises.

Pourquoi les marées sont importantes pour la voile?

Les deux questions génériques pour la navigation à voile est d’identifier quel sera le niveau d’eau à un moment donné et sa question inverse, à savoir à quelle heure aura lieu un niveau d’eau donné. Toutes les questions de navigation liées à l’ancrage, aux moments de passages, de même que les questions usuelles sur les courants et les niveaux d’eau reviennent à ces deux questions. Une approche moderne de navigation consiste à employer une application de navigation, mais la prudence demanderait d’avoir également une approche alternative qu’on puisse faire à bord.

L’approche de la veille école pour répondre aux deux questions génériques est de se fier aux tables des marées. Pour un (demi) cycle de marée, on produit une interpolation manuelle, permettant ensuite de répondre aux questions. Pour cette raison, apporter une copie de la table des marées à bord est l’approche standard pour d’éventuelles adaptations à un plan de navigation.

Les techniques d’interpolation varient selon les écoles de navigation. J’ai déjà écrit sur la méthode du RYA. Leur technique s’appuie sur une approche graphique (non-paramétrique) qui représente un cycle de marée typique. La technique du RYA fonctionne très bien quand les marées ont structurellement la même forme (comme dans le Solent). L’approche a l’avantage d’incorporer les déviations réelles à une courbe sinusoïdale dans les prévisions. Par contre, si les marées ne sont pas structurellement stables, comme c’est le cas des marées mixtes, cette approche devient moins attrayante car elle demande plusieurs graphiques pour s’adapter aux circonstances.

Dans un billet précédent, j’avais souligné que les deux méthodes d’interpolation enseignées par Voile Canada produisent sensiblement les mêmes résultats en dépit des apparences différentes. La fameuse « règle des sixièmes » et les tableaux 5 et 5A de Pêches et Océans sont des approximations d’une courbe sinusoïdale. Les deux produiront des interpolations similaires, avec des écarts de plus ou moins cinq centimètres. J’avais conclu en évoquant la possibilité d’apporter un modèle de prévision de marée en guise d’alternative aux volumineuses tables des marées.

L’exercice que je suis sur le point de présenter m’a convaincu d’apporter à bord 37 nombres. Ce sont les seuls intrants nécessaires pour faire des prévisions de marées à n’importe quel endroit au monde… pour autant que des données sur les niveaux d’eau existent. De plus, si les prévisions requises sont pour le futur immédiat, il est possible de prédire le niveau des marées avec une précision accrue. Finalement, la technique d’estimation fonctionne pour n’importe quel moment de la journée, éliminant de fait le besoin d’interpolation.

Mais aussi attrayante qu’elle puisse paraître, la technique ne m’a pas convaincue d’abandonner les méthodes « usuelles » de calculs des marées. D’abord, l’approche est exigeante sur le plan technique: il faut écrire du code informatique. C’est beaucoup plus simple d’utiliser le code des autres à travers une application. Et pour une estimation rapide, entre deux tâches sur le pont, la règle des sixièmes est difficile à battre.

J’ai appris cette technique de prévision en me fiant sur quatre références, soit le document de référence de la National Oceanic and Atmospheric Administration (NOAA; Parker, 2007), un article de recherche par Foremann, Cherniwawsky et Ballantine (2009) sur les prévisions de marées, un chapitre de référence sur les prévisions de marées (Foremann, Crawford et Marsden, 1995) et la page Wikipedia (!) sur l’analyse spectrale par moindre carrés (Wikipeia, s.d.). Le chapitre 3 de Parker est particulièrement utile pour faire le pont entre la physique et les techniques statistiques de prévision.

La suite est pleine de maths

Je ne prétends pas que cette méthode est pour tout le monde. Il faut une compréhension des mathématiques et des statistiques qui est soit de premier cycle universitaire, soit de cycles supérieurs universitaires. Il faut comprendre les séries de Fourier et les régressions linéaires. Les séries de Fourier sont des outils mathématiques particulièrement adaptés pour représenter des phénomènes cycliques (comme les marées). Les régressions linéaires sont des outils statistiques qui permettent d’identifier quels modèles mathématiques reproduisent le mieux les données réelles.

C’est donc possible de lire ce texte pour le plaisir, en laissant les prévisions aux matheux. C’est aussi possible de faire les prévisions vous-même. Tout ce qui est présenté ci-dessous peut se faire avec des logiciels gratuits. Le choix vous appartient, mais vous êtes avertis: la suite est pleine de maths!

Modélisation des marées 101

Les phénomènes périodiques peuvent se modéliser en fonction du temps t à partir d’une somme de cosinus. C’est un cas particulier de série de Fourier. Une représentation canonique d’un modèle ressemble à:

h(t) = \beta_0 + \sum_{i=1}^n H_i\cos(\alpha_i t+k_i),

où h(t) est le phénomène d’intérêt (p. ex. la hauteur de la marée) et où \beta_0, H_i, \alpha_i et k_i, i=1...n sont des paramètres inconnus. Choisir ces paramètres de manière appropriée est la clé pour reproduire le phénomène cyclique. Dans certains cas, cette décomposition en série de cosinus n’est rien d’autre qu’un exercice de modélisation. Augmenter le nombre de cosinus augmente la précision prévisionnelle (voir le vidéo ci-dessous) et réduit l’imprécision de prévision.

Dans le cas des marées, les cosinus ont cependant un sens physique. Ils représentent un effet spécifique découlant de l’orbite de la lune et du soleil sur l’eau. Si l’orbite de la terre et de la lune étaient dans le même plan, il ne faudrait que deux cosinus pour prédire les marées (n serait égal à deux). Mais ces orbites sont elliptiques et ne sont pas dans le même plan, si bien qu’elles génèrent des « non-circularités » qu’on doit représenter par des cosinus additionnels. Comme le titre de ce billet le suggère, 37 cosinus font un bon travail pour fins de modélisation (n=37), mais c’est possible d’en ajouter d’autres.

Historiquement, la prévision des marées avait beaucoup à voir avec l’identification des inconnus H_i, \alpha_i et k_i à partir des principes de physique. Les inconnus \alpha_i sont maintenant bien connus et sont désormais des « constantes » liées à la fréquence des orbites de la lune et de la terre (ou de leur combinaison). Ils sont mesurés en degrés par heure et réfèrent à la vitesse à laquelle des phénomènes se répètent. Les autres composantes dépendent de phénomènes locaux. Ils étaient originalement décomposés en phénomènes généraux et locaux, puis ajustés manuellement pour améliorer les prévisions.

Il y a cependant un très bon truc pour éviter de trouver ces constantes à partir d’une modélisation physique. Avec l’usage d’une identité trigonométrique, l’expression h(t) peut être transformée en une expression linéaire des inconnues. Cette expression linéaire peut ensuite être estimée par régression linéaire, ce qui permet d’identifier les coefficients à partir des données. Bref, au lieu de se casser la tête avec un modèle de physique, on laisse les données de marée révéler la valeur des inconnus. L’approche est toujours utilisée aujourd’hui et s’appelle (traduction libre) un modèle harmonique par moindres carrés.

Comment estimer les coefficients de marée

On se rappelle que le cosinus de deux angles peut être décomposé par une somme de produits de sinus et de cosinus:

\cos(a+b) = \cos(a)\cos(b)+\sin(a)\sin(b).

En utilisant cette expression pour notre modèle de marée h(t), on peut écrire:

\cos(\alpha_i t +k_i) = \cos(\alpha_i t)\cos(k_i) + \sin(\alpha_i t)\sin(k_i).

Le modèle devient alors:

h(t) = \beta_0 + \sum_{i=1}^n\left(\cos(\alpha_i t)\underbrace{H_i\cos(k_i)}_{\equiv\beta_{1i}}+ \sin(\alpha_i t)\underbrace{H_i\sin(k_i)}_{\equiv\beta_{2i}}\right).

On se souvient que H_i et k_i sont deux inconnus, si bien que chaque produit H_i\cos(k_i), H_i\sin(k_i), peut être vu comme deux nouvelles inconnues. En les notant \beta, on obtient:

h(t) = \beta_0 + \sum_{i=1}^n\left(\cos(\alpha_i t)\beta_{1i} + \sin(\alpha_i t)\beta_{2i}\right).

Si les coefficients \alpha_i sont connus, les expressions de sinus et de cosinus sont des nombres connus, si bien que la fonction h(t) est une combinaison linéaire des inconnues \beta. On peut les identifier par régression linéaire.

Pour des échantillons de grande taille, trouver une méthode d’estimation « assez rapide » est d’importance pratique. Certaines techniques d’estimation fonctionnent bien en théorie ou sur de petits échantillons, mais prennent tellement de temps sur des grands échantillons qu’elles deviennent inutilisables. Une borne supérieure pratique est que le temps requis de l’exécution d’une méthode d’estimation devrait être inférieur à la fenêtre d’estimation qu’elle fournit. En d’autres termes, un modèle qui fait des prévisions pour les trois prochains jours devrait prendre un temps d’estimation… de moins de trois jours!

À cet égard, être capable de convertir une somme de cosinus en expression linéaire est point crucial d’usage pratique. Les régressions linéaires s’estiment rapidement, elles sont robustes et elles se moquent complètement d’éventuelles données manquantes. Sur un ordinateur personnel, on peut les estimer en une fraction de seconde. Ce n’est pas le cas de techniques alternatives servant à estimer des cosinus. Cette conversion est possible parce qu’on connaît le mouvement orbital de la terre et de la lune, ce qui permet de calculer les \alpha_i. Comme c’est le cas de la plupart des technologies applicables, c’est une combinaison de connaissances théoriques et de méthodes empiriques qui fournit une approche fonctionnelle.

Quels sont les 37 nombres?

Tiré de Parker (2007, pp 88-89):

Ce n’est pas l’ambition de ce texte que d’expliquer l’origine de chacun de ces nombres. On peut lire Parker (2007) ou Foremann et al. (1995) pour une explication détaillée. Pour les besoins de prévisions, ces 37 nombres sont à retenir.

Deux exemples

J’utilise la méthode de prévision harmonique par moindres carrés pour faire des prévisions de marées à deux endroits au Canada: Rimouski (Québec) et Point Atkison (Colombie-Britannique). J’ai choisi ces deux endroits parce que faire un calcul manuel de marée à Rimouski est un classique d’examen de Voile Canada et parce que Pêches et Océans Canada m’a gracieusement fourni les données historiques de Point Atkinson (merci!).

Pêches et Océans Canada partage publiquement ses donnée, bien que parfois avec un historique limité et seulement pour certaines stations de mesures. Ils ont une API publique (ici) à partir duquel on peut télécharger les données. J’ai téléchargé cinq années de données pour Rimouski. Bien que l’échantillon soit de petite taille, il est à jour, au sens où la dernière observation est à la même date que la journée de rédaction de ce billet.

L’échantillon pour Point Atkinson possède 25 années de données (avec des donnée manquantes), mais la dernière observation est en 2019. Les données sont assujetties à la déclaration écrite (en anglais) à l’Annexe 1 et le code Python pour la récupération de données est à l’Annexe 2. Un résumé des deux jeux de données est disponible ci-dessous.

Rimouski, QC

Les dernières 500 heures de données sont illustrées dans la Figure ci-dessous.

L’essentiel du code informatique pour estimer le modèle est détaillé ci-dessous. La ligne 35 définit les 37 nombres discutés ci-dessus. Les lignes 39 à 45 définissent les cosinus et sinus du modèle de marée. Les lignes 47 et 48 font la régression linéaire, alors que la ligne 49 illustre les résultats de la régression.

Le modèle résultant de cette estimation est affiché (partiellement) ci-dessous. L’image illustre deux encadrés. Le premier en haut à gauche montre un R-carré de 0.97, montrant que le modèle reproduit 97% de la variance contenue dans les données. Le tableau au bas de l’image montre les coefficients estimés du modèle (les \beta_i) de même que les informations requises pour faire les tests statistiques de significativité. Certaines lignes de ce tableau montrent une valeur p (encadré rouge du bas, au centre) supérieure à 5%, ce qui suggère que les cosinus/sinus associés apportent davantage de bruit que de signal au modèle de prévision (les statistique t et les intervalles de confiance fournissent la même information, mais présentée de manière différente). Bref, un coefficient qui apporte plus de bruit que d’aide à la prévision devrait être retiré pour réduire la marge d’erreur prévisionnelle.

Quels coefficients devraient être conjointement retirés se calcule à partir d’une statistique de Fischer spécifique, parce que les valeurs p sont interdépendantes. Cela dit, les régressions linéaires sont tellement rapides à estimer qu’il est plus commode de faire (et d’automatiser) le retrait d’un coefficient, la ré-estimation du modèle résultant, puis d’évaluer quels autres coefficients sont à retirer.

Ce processus itératif peut-être répété jusqu’à ce que tous les coefficients estimés résultants soient statistiquement différents de zéro. Le modèle résultant de cette méthodologie (qui n’est pas illustré) affiche toujours un R-carré de 0.97, montrant que les variables retirées apportent peu en pouvoir explicatif. L’erreur quadratique moyenne est de 0.027, c’est-à-dire un peu moins de trois centimètres d’erreur de prévision.

La prévision illustrée ci-dessus montre que le modèle a plus de difficultés avec les marées mortes (la ligne bleu et la ligne jaune diffèrent). Les prévisions sont définitivement très proches de celles de Pêches et Océans Canada et un examen rapide de ce que fait marees.gc.ca suggère également que leur modèle est moins performant pendant les marées mortes (détails plus bas). Le tableau ci-dessous sera rempli au fur et à mesure que les marées se réalisent.

Point Atkinson, BC

Point Atkinson a des marées mixtes. Les dernières 500 heures sont présentées ci-dessous. Je présente également le modèle de prévision découlant de la même méthodologie (code en annexe). L’absence de commentaire de cette section est une éloge implicite quant à la généralité de l’approche harmonique, car il n’y a aucun changement à la méthodologie.

Peut-on rendre l’approche plus performante?

Une seule équation estimée à partir d’une technique statistique simple peut fournir des prévisions de marées pour des périodes de temps arbitrairement longues. Ce progrès technologique impressionnant et simple est le résultat du travail scientifique cumulé de plusieurs personnes sur une période de plusieurs siècles. L’approche est extrêmement utile et donne une bonne idée de la manière dont on fait les prévisions contenues dans les tables des marées.

Pourtant, un examen des différences entre les valeurs réelles et les valeurs prédites suggère qu’il y a parfois des différences importantes, particulièrement pendant les marées mortes. Dans certains cas, la différence va jusqu’à 20 centimètres (8 pouces) lors d’un extrema. Comment peut-on améliorer la performance prévisionnelle?

Une amélioration facile consiste à utiliser les erreurs de prévision du modèle harmonique pour tenter de prédire les erreurs de prévisions futures. On se rappelle que les marées sont le fruit des forces gravitationnelles, mais aussi des effets locaux. L’inclusion de méthodes de prévisions d’erreurs passées peut donc se justifier par le fait que le modèle original n’incorpore pas les conditions locales affectant les marées. Mathématiquement, on peut inclure les erreurs de prévision en introduisant une composante auto-régressive au modèle:

h(t) =\beta_0+ \sum_{i=1}^n\left(\cos(\alpha_it)\beta_{1i} + \sin(\alpha_it)\beta_{2i} \right) + \sum_{l=1}^L\rho_l h(t-l).

La dernière sommation inclut des valeurs passées du niveau d’eau, jusqu’à concurrence de L retards et leur effet sur les niveaux présents sont traduits par les coefficients \rho_l qui sont à estimer.

Cette approche peut drastiquement améliorer la performance prévisionnelle, mais a un prix. Puisqu’il faut les données réalisées du passé, il faut que ce passé se soit réalisé pour être en mesure de faire la prévision. En d’autres termes, si la dernière heure est requise pour fin de prévisions, on ne peut pas faire de prévisions de plus qu’une heure d’avance. Conséquemment, ce genre de modèle vient avec un arbitrage à faire entre le besoin de prédire à des périodes futures qui sont éloignées vis-à-vis une précision accrue.

Le point de départ de l’analyse auto-régressive consiste à évaluer quelles valeurs passées sont les plus corrélées avec les valeurs présentes. On peut faire cette analyse à l’aide d’un diagnostic d’autocorellation partielle sur les erreurs de prévision. Les deux lignes de code ci-dessous se chargent de produire le graphique d’intérêt (sur les données de Rimouski):

La ligne 63 calcule les erreurs de prévision alors que la ligne 64 produit le graphique d’autocorrélation partielle (ci-dessous). Bien que ce ne soit pas apparent, la marge d’erreur de 95% déterminant le seuil de rejet de l’hypothèse nulle (aucune corrélation) est collé sur l’axe horizontal du graphique.

Ces résultats signifient qu’il y a énormément d’information qui est dans les erreurs de prévision et qui n’est pas exploitée. Leur usage permettrait donc d’améliorer les prévisions. Ce qui devrait sauter aux yeux du graphique ci-dessous est la capacité prédictive des trois dernières heures de marées, représentées par des pics élevés sur la gauche du graphique. Ces derniers sont assurément les meilleures améliorations pour fins de prédictions… mais condamnent également le modèle à ne faire que des prévisions à l’heure.

L’inclusion de niveaux d’eau passé dans le modèle transforme la technique d’estimation statistique. On passe des modèles de régression linéaires aux modèles de régressions arimax. Normalement, ces modèles sont estimés par méthode de maximum de vraisemblance, mais un peu de théorie statistique (que je ne détaille pas ici) nous informe que l’estimation par moindre carrés sera également consistante si l’échantillon est d’assez grande taille.

J’emploie cette deuxième approche dans l’extrait de code ci-dessous. Les lignes 69 à 74 créent les variables décalées alors que les lignes 77 à 80 refont l’estimation du modèle.

La prévision découlant de ce modèle est présentée dans la Figure ci-dessous. Il n’y a maintenant plus de différence pratique entre les données réelles et les prévisions. Si le but est de reproduire les données, ou encore de faire des prévisions une heure à l’avance, ce modèle est très difficile à battre. Il n’est cependant d’aucune aide pour des prévisions à plus longue échéance.

Conclusion

Le fait que 37 nombres permettent de construire un modèle performant de prévision des marées est impressionnant. Étant donné que ce modèle est central aux documents de la NOAA, il donne une assez bonne impression de la manière dont les tables des marées sont construites. Je ne serais pas surpris si ce genre de modèle était imbriqué dans un ordinateur de bord.

C’est peut-être une singularité des deux jeux de données retenus pour fins d’illustration, mais ce modèle semble être moins performant pendant les marées mortes. Un bon complément prévisionnel à ce modèle serait de faire des prévisions à l’heure lors de marées mortes (comme dans la dernière section).

Les modèles de court-terme de Pêches et Océans Canada permettent de faire une prévision jusqu’à une semaine d’avance. Leur échantillonnage est à la minute. Ils présentent un modèle harmonique pour les prévisions de long-terme, mais leur capacité à faire une prévision sur une semaine n’est pas quelque chose que j’ai élucidé… et me fera peut-être revenir sur la question de l’estimation des marées.

Est-ce qu’un modèle harmonique peut-être un substitut à une table des marées? Peut-être. Un petit modèle de cette forme comprend beaucoup plus d’information qu’une table de prévision des marées pour la même région… si les données sont disponibles. Le fait qu’il ne se base que sur 37 constantes le rend également élégant. Je peux concevoir l’usage d’un tel modèle dans une région du monde où il y a des données de niveau d’eau, mais aucune prévision de disponible. Pour un ajustement rapide aux plans de navigation, je me fierais cependant sur une application de navigation avec une réputation établie. C’est plus facile d’usage.

Remerciements

Je tiens à remercier Mike Foreman pour m’avoir aidé à identifier des références sur les modèles de prévisions de marées. Toutes les erreurs demeurent les miennes.

Références

3Blue1Brown (2018). But What is a Fourier Transform? From Heat Flow to Drawing With Circles, vidéo YouTube récupéré en ligne en novembre 2023 à this address.

Admiralty Maritime Data Solutions (n.d.). Calshot Castle, England [Tide Data], dynamic webpage retrieved in Novembre 2023 at this address.

Pêches et Océans Canada (n.d.). Rimouski (02985) [Tidal Station Information], page web dynamique récupérée en ligne en novembre 2023 à this address.

_________________________ (n.d.). API for the Integrated Water Level System, page web récupérée en ligne en novembre 2023 à this address.

Foremann, M. G. G., Cherniawsky, J. Y. et V. A. Ballantyne (2009). Versatile Harmonic Tidal Analysis: Improvements and Applications, Journal of Atmospheric and Oceanic Technology V(26), pp 806-817.

Foremann, M. G. G., Crawford, W. R. et R. F. Marsden (1995). De-Tiding: Theory and Practice, Chapter 10 of Coastal and Estuarine Studies V(47), pp 203-239.

Jeandusud.com (2023). RYA YachtMaster et Programme de quillard de Voile Canada: 11 différences, text availlable at this address.

_____________ (2023). La règle des sixièmes et des vingtièmes, text availlable at this address.

Parker, B. B. (2007). Tidal Analysis and Prediction, NOAA Special Publication NOS CO-OPS 3, récupérée en ligne en novembre 2023 à at this address.

Waterlust (2022). How the Tides REALLY Work, vidéo YouTube récupéré en ligne en novembre 2023 à this address.

Wikipedia (n.d.). Séries de Fourier, page web récupérée en ligne en novembre 2023 à this address.

________ (n.d.). Least-Square Spectral Analysis, page web récupérée en ligne en novembre 2023 à this address.

________ (n.d.). Régressions linéaires, page web récupérée en ligne en novembre 2023 à this address.

________ (n.d.). Liste des identités trigonométriques, page web récupérée en ligne en novembre 2023 à this address.

________ (n.d.). Partial Auto-correlation Function, page web récupérée en ligne en novembre 2023 à this address.

Annexe 1: Notice du Gouvernement du Canada

This product is not to be used for navigation.

This product was made by or for Jeandusud.com and contains intellectual property (Data) of the Canadian Hydrographic Service of the Department of Fisheries and Oceans.

The copyrights in the Data are and remain the property of His Majesty the King in Right of Canada and shall not be sold, licensed, leased, assigned or given to a third party.

The incorporation of the Data in this product does not constitute an endorsement or an approval of this product by the Canadian Hydrographic Service, the Department of Fisheries and Oceans or His Majesty the King in Right of Canada.

Annexe 2: code informatique

Récupération de données (Rimouski)

from urllib.request import urlopen 
import json, time
import pandas as pd

#Rimouski
api_url_1 = 'https://api-iwls.dfo-mpo.gc.ca/api/v1/stations/5cebf1e03d0f4a073c4bbd92/data?time-series-code=wlo'
api_url_2 = '&resolution=SIXTY_MINUTES'
for some_year in range(2019, 2024):
	if some_year == 2019: 
		month_range = range(8, 13)
	elif some_year == 2023: 
		month_range = range(1, 12)
	else:
		month_range = range(1, 12)
	for some_month in month_range: 
		if some_year < 2023 or (some_year == 2023 and some_month <= 10):
			if some_month != 12: 
				complete_url = api_url_1 + "&from=" + str(some_year) + "-" + str("%02d" % (some_month,)) + "-" +  "27T00%3A00%3A00Z" + "&to=" + str(some_year) + "-" + str("%02d" % (some_month+1,)) + "-" + "27T00%3A00%3A00Z" + api_url_2
			else:
				complete_url = api_url_1 + "&from=" + str(some_year) + "-" + "12" + "-" +  "27T00%3A00%3A00Z" + "&to=" + str(some_year+1) + "-" + "01" + "-" + "27T00%3A00%3A00Z" + api_url_2
		response = urlopen(complete_url) 
		tidal_data = json.loads(response.read())
		if some_year == 2019 and some_month == 8:
			data = pd.DataFrame(tidal_data)
		else:
			data_temp = pd.DataFrame(tidal_data)
			data = data.append(data_temp, ignore_index=True)
		time.sleep(1)
data['eventDate'] = pd.to_datetime(data['eventDate'])
data.to_csv("rimouski.csv")

Analyse (Rimouski)

import pandas as pd
import numpy as np
import math
from scipy.signal import argrelextrema
from datetime import datetime
from datetime import timedelta
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import statsmodels.api as sm

data = pd.read_csv('rimouski.csv')
data['eventDate'] = pd.to_datetime(data['eventDate'])
data['t'] = data['eventDate'] - data['eventDate'][0]
data['seconds'] = data['t'].apply(timedelta.total_seconds)
data['hours'] = data['seconds']/3600
original_size = data.shape[0]

fig, axs = plt.subplots(figsize=(20, 4))
axs.plot(data.loc[data.index > data.shape[0]-501, 'eventDate'], data.loc[data.index > data.shape[0]-501, 'value'])
plt.xlabel("Date")
plt.ylabel("Level (m)")
plt.show()

#Expanding dataset for forecast
exp = pd.DataFrame(pd.date_range(data['eventDate'][data.shape[0]-1], periods=101, freq='h'))
exp.columns = ['eventDate']
exp = exp[1:]
exp["hours"] = 1
exp["hours"][0] = data['hours'][data.shape[0]-1]+1
exp["hours"] = exp["hours"].cumsum()
exp["value"] = np.nan

data = data.append(exp, ignore_index=True)


period_numbers = [28.9841042, 57.9682084, 86.9523127, 115.9364169, 15.0410686, 90, 13.9430356, 16.1391017, 42.9271398, 44.0251729, 12.8542862, 30, 31.0158958, 58.9841042, 60, 14.4920521, 43.4761563, 28.4397295, 13.3986609, 15.5854433, 57.4238337, 29.5284789, 27.9682084, 0.5443747, 14.9589314, 30.0821373, 1.0158958, 1.0980331, 28.5125831, 29.4556253, 27.8953548, 13.4715145, 0.0410686, 0.0821373, 15, 30.0410667, 29.9588333]
data['b0'] = 1
variables = ['b0']

for i in range(0, len(period_numbers)):
	cf = lambda t, j=i : math.cos(2*math.pi/360*t*period_numbers[j])
	sf = lambda t, j=i : math.sin(2*math.pi/360*t*period_numbers[j])
	data['b1' + str(i)] = data['hours'].apply(cf)
	data['b2' + str(i)] = data['hours'].apply(sf)
	variables.append('b1'+str(i))
	variables.append('b2'+str(i))

for some_var in ['b13','b23', 'b15', 'b25', 'b18', 'b19' , 'b210', 'b214', 'b112', 'b116', 'b216', 'b221', 'b134','b235']:
	variables.remove(some_var)

model = sm.OLS(data['value'][0:original_size], data[variables][0:original_size])
results = model.fit()
print(results.summary())
data['forecast'] = results.predict(data[variables])

fig, axs = plt.subplots(figsize=(20, 4))
axs.plot(data.loc[data.index > data.shape[0]-501, 'eventDate'], data.loc[data.index > data.shape[0]-501, 'value'])
axs.plot(data.loc[data.index > data.shape[0]-501, 'eventDate'], data.loc[data.index > data.shape[0]-501, 'forecast'])
plt.xlabel("Date")
plt.ylabel("Level (m)")
plt.legend(['Observed', 'Predicted'])
plt.show()

data['errors'] = data['forecast'] - data['value']
w = np.linspace(0.01, 10, 10000)
f = scipy.signal.lombscargle(data['hours'][0:original_size], data['errors'][0:original_size], w)
plt.plot(w, f)


data['value_L1'] = data['value'].shift(1)
data['value_L2'] = data['value'].shift(2)
data['value_L3'] = data['value'].shift(3)
variables.append('value_L1')
variables.append('value_L2')
variables.append('value_L3')


model = sm.OLS(data['value'][28:original_size], data[variables][28:original_size])
results = model.fit()
print(results.summary())
data['forecast2'] = results.predict(data[variables])

fig, axs = plt.subplots(figsize=(20, 4))
axs.plot(data.loc[data.index > data.shape[0]-501, 'eventDate'], data.loc[data.index > data.shape[0]-501, 'value'])
axs.plot(data.loc[data.index > data.shape[0]-501, 'eventDate'], data.loc[data.index > data.shape[0]-501, 'forecast2'])
plt.xlabel("Date")
plt.ylabel("Level (m)")
plt.legend(['Observed', 'Predicted'])
plt.show()

w = np.linspace(0.01, 10, 10000)
f = scipy.signal.lombscargle(data['hours'][0:original_size], data['errors'][0:original_size], w)
plt.plot(w, f)
#big spike at 6,24

Analyse (Point Atkinson)

import pandas as pd
import numpy as np
import math
from scipy.signal import argrelextrema
from datetime import datetime
from datetime import timedelta
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import statsmodels.api as sm

data = pd.read_csv('7795-1995-2019_hrly-slev.csv', sep=";")
columns = data.columns.tolist()
columns[2] = "day"
columns[3] = "hour"
columns[4] = "minute"
data.columns = columns
data["t"] = pd.to_datetime(data[columns[0:4]])
data['t'] = pd.to_datetime(data['t'])
data['delta'] = data['t'] - data['t'][0]
data['seconds'] = data['delta'].apply(timedelta.total_seconds)
data['hours'] = data['seconds']/3600
original_size = data.shape[0]

fig, axs = plt.subplots(figsize=(20, 4))
axs.plot(data.loc[data.index > data.shape[0]-501, 't'], data.loc[data.index > data.shape[0]-501, 'level'])
plt.xlabel("Date")
plt.ylabel("Level (ft)")
plt.show()

#Expanding dataset for forecast
exp = pd.DataFrame(pd.date_range(data['t'][data.shape[0]-1], periods=100, freq='h'))
exp.columns = ['t']
exp["hours"] = 1
exp["hours"][0] = data['hours'][data.shape[0]-1]+1
exp["hours"] = exp["hours"].cumsum()
exp["level"] = np.nan

data = data.append(exp, ignore_index=True)


period_numbers = [28.9841042, 57.9682084, 86.9523127, 115.9364169, 15.0410686, 90, 13.9430356, 16.1391017, 42.9271398, 44.0251729, 12.8542862, 30, 31.0158958, 58.9841042, 60, 14.4920521, 43.4761563, 28.4397295, 13.3986609, 15.5854433, 57.4238337, 29.5284789, 27.9682084, 0.5443747, 14.9589314, 30.0821373, 1.0158958, 1.0980331, 28.5125831, 29.4556253, 27.8953548, 13.4715145, 0.0410686, 0.0821373, 15, 30.0410667, 29.9588333]
data['b0'] = 1
variables = ['b0']

for i in range(0, len(period_numbers)):
	cf = lambda t, j=i : math.cos(2*math.pi/360*t*period_numbers[j])
	sf = lambda t, j=i : math.sin(2*math.pi/360*t*period_numbers[j])
	data['b1' + str(i)] = data['hours'].apply(cf)
	data['b2' + str(i)] = data['hours'].apply(sf)
	variables.append('b1'+str(i))
	variables.append('b2'+str(i))

model = sm.OLS(data['level'][0:original_size], data[variables][0:original_size])
results = model.fit()
print(results.summary())

for some_var in ['b13','b23', 'b15', 'b25', 'b18', 'b29' , 'b114', 'b214', 'b116', 'b216', 'b131', 'b135','b235']:
	variables.remove(some_var)
data['forecast'] = results.predict(data[variables])

fig, axs = plt.subplots(figsize=(20, 4))
axs.plot(data.loc[data.index > data.shape[0]-501, 't'], data.loc[data.index > data.shape[0]-501, 'level'])
axs.plot(data.loc[data.index > data.shape[0]-501, 't'], data.loc[data.index > data.shape[0]-501, 'forecast'])
plt.xlabel("Date")
plt.ylabel("Level (m)")
plt.legend(['Observed', 'Predicted'])
plt.show()