Terschelling (Netherlands)¶
import warnings
warnings.filterwarnings("ignore")
# os
import os
import os.path as op
import sys
# arrays
import math
import numpy as np
import pandas as pd
from scipy import signal as sg
# plot
import matplotlib.pyplot as plt
sys.path.insert(0, os.path.join(os.getcwd() , '..', '..', '..'))
# dependencies
if(os.path.isdir('waves-main')): #thebe
os.chdir('waves-main')
from lib.eta_spec import *
# path to data
if(os.path.isdir('data')):
p_data = op.abspath(op.join(os.getcwd(), 'data', 'Terscheling')) # thebe
else:
p_data = op.abspath(op.join(os.getcwd(),'..', '..', '..', 'data', 'Terscheling')) # notebook
fop04930 = np.loadtxt(op.join(p_data, "fop04930.txt"))
fop05235 = np.loadtxt(op.join(p_data,"fop05235.txt"))
fra05092 = np.loadtxt(op.join(p_data,"fra05092.txt"))
fre05400 = np.loadtxt(op.join(p_data,"fre05400.txt"))
fri05276 = np.loadtxt(op.join(p_data,"fri05276.txt"))
fri05658 = np.loadtxt(op.join(p_data,"fri05658.txt"))
# Select deploy
datos = fop04930
Field data
Sampling frequency: 2 seconds
A.1 Observations¶
# The free surface elevation can be seen as a Gaussian process
plt.figure(figsize=(13,4))
plt.plot(datos, color='b')
plt.xlabel('Time (2s sampling interval)')
plt.ylabel('$\eta$ (m)')
plt.xlim([0, len(datos)])
plt.show()
print("Mean: ", np.mean(datos))
print("Standar Deviation: ", np.var(datos))
Mean: 2.197265624999758e-07
Standar Deviation: 0.014529366594190001
A.2. Spectral Analysis¶
# Calculate spectra help(sg.welch)
f, E = sg.welch(datos, fs = 2, nfft = len(datos))
# Plot wave spectrum
plt.figure(figsize=(5,5))
plt.plot(f, E, c='b')
plt.xlabel('frequency (Hz)')
plt.ylabel('Densidad espectral (/Hz)')
plt.ylabel('PSD (m²/Hz)')
plt.show()
Wave spectrum parameters¶
\[ m_{n} = \int_{0}^{\infty}S(\omega)d\omega ; n=0,1,2... \]
\(m_{0}:\) zeroth-order moment of the wave spectrum
\(m_{1}:\) first-order moment of the wave spectrum
# Momento calculado con el método calc_momento
print("Zeroth-order moment: " + str(moment(0, f, E)))
print("First-order moment: " + str(moment(1, f, E)))
Zeroth-order moment: 0.013835893119162548
First-order moment: 0.0025636608788534593
\[ H_{m0} = 4.004 \sqrt{m_{0}} \]
# Wave height of the zeroth order moment
print("Wave heigh: " + str(4.004 * math.sqrt(moment(0, f, E))))
Wave heigh: 0.47097479747933846
\[ \overline{\omega}=\frac{m_{1}}{m_{0}};\; \overline{T}=2\pi\frac{m_{0}}{m_{1}};\; \overline{T_{z}}=2\pi\frac{m_{0}}{m_{2}} \]
# Mean frequency
print("Mean frequency : " + str(moment(1, f, E)/moment(0, f, E)))
# Wave period associated to mean frequency
print("Wave period associated to mean frequency: " + str(moment(0, f, E)/moment(1, f, E)))
# Mean zero-crossing period
print("Mean zero-crossing period: " + str(np.sqrt(moment(0, f, E)/moment(2, f, E))))
Mean frequency : 0.18529059575509582
Wave period associated to mean frequency: 5.396927976429685
Mean zero-crossing period: 5.170505554148997
\[\overline{T_{e}}=T_{-1}=2\pi\frac{m_{-1}}{m_{0}}; \;
v^{2}=\frac{m_{0}m_{2}}{m_{1}^{2}}-1
\]
# Mean energy period
print("Mean energy period: " + str(moment(-1, f, E)/moment(0, f, E)))
# Spectral width of Longuet-Higgins (1957)
# When de energy is concentrated in one single frequency, v = 0
print("Spectral width: " + str((moment(0, f, E) * moment(2, f, E)/moment(1, f, E)**2)-1))
Mean energy period: 5.947513682876696
Spectral width: 0.08949998015393712
JONSWAP Goda (1985)¶
\(S(f)=\alpha\cdot H_{s}^{2} \cdot T_{p}^{-4} \cdot f^{-5} \cdot e^{-1.25 \cdot (T_{p} \cdot f)^{-4}} \cdot \gamma^{e^{-(T_{p} \cdot f - 1)^{2}/(2 \cdot \sigma^{2})}}\)
\(\alpha \approxeq \frac{0.0624}{0.230+0.0336 \cdot \gamma - 0.185 \cdot (1.9 + \gamma)^{-1}}\)
\( \sigma=\begin{cases} \sigma_{a}; f \leq f_{p} \\ \sigma_{b}; f \geq f_{p} \end{cases} \)
\(\gamma = 1 to 7 (mean 3.3), \sigma_{a}\approxeq0.07, \sigma_{b}\approxeq 0.09\)
print('Peak frequency: ' + str(f[np.argmax(E)]))
print('Peak period: ' + str(1/f[np.argmax(E)]))
Peak frequency: 0.125
Peak period: 8.0
gamma, EJon = assess_jonwsap(f, E)
print('Best gamma JONSWAP fit: ' + str(gamma))
Best gamma JONSWAP fit: 2
# Plot measured wave spectrum and theorecial JONSWAP shapes
plt.figure(figsize=(5,5))
plt.plot(f, E, c='b', label='Measured')
plt.plot(f[:-1], EJon, c='r', label='JONSWAP gamma = {0}'.format(gamma))
plt.xlabel('frequency (Hz)')
plt.ylabel('Densidad espectral (/Hz)')
plt.ylabel('PSD (m²/Hz)')
plt.legend()
plt.show()
A.3 Short-term statistics¶
Definition of ‘waves’ in a time record of the surface elevation with upward zero-crossings
fs = 1/2
T, H = upcrossing(datos, fs)
fig, axs = plt.subplots(2, 1, figsize=(15,5), sharex=True)
axs[0].vlines(range(len(H)), np.full(len(H), 0), H, color='orangered')
axs[1].vlines(range(len(T)), np.full(len(T), 0), T, color='slateblue')
axs[0].scatter(range(len(H)), H, color='orangered')
axs[1].scatter(range(len(T)), T, color='slateblue')
axs[0].set_ylabel('H (m)')
axs[1].set_ylabel('T (s)')
axs[1].set_xlabel('sequential waves in a sea-state (1-hour)')
plt.xlim(0, 100)
plt.show()
Statistical parameters¶
mean wave heigh \(\overline H\)
\(\overline H = \frac{1}{N} \sum_{i=1}^{N} H_{i}\)
where i is the sequence number (in time) of the wave in the record
root-mean-square wave height \(H_{rms}\)
\(H_{rms}=(\frac{1}{N} \sum_{i=1}^{N} H_{i}^{2})^{1/2}\)
significant wave heigh \(H_{1/3}\)
\(H_{1/3}=\frac{1}{N/3} \sum_{j=1}^{N/3} H_{j}\)
where j is the rank number of the wave, based on wave-heigh
mean of the highest one-tenth of waves \(H_{1/10}\)
\(H_{1/10}=\frac{1}{N/10} \sum_{j=1}^{N/10} H_{j}\)
mean zero-crossing wave period \(\overline T_{0}\)
\(\overline T_{0}=\frac{1}{N} \sum_{i=1}^{N} T_{0,i}\)
significant wave period \(T_{1/3}\)
\(T_{1/3}=\frac{1}{N/3} \sum_{j=1}^{N/3} T_{0,j}\)
# mean wave heigh
print("mean wave heigh: " + str(np.mean(H)))
# root-mean-square wave height
print("root-mean-square wave height: " + str(rmsV(H)))
# significant wave heigh
print("significant wave heigh: " + str(highestN_stats(H, 3)))
# mean of the highest one-tenth of waves
print("mean of the highest one-tenth of waves: " + str(highestN_stats(H, 10)))
# maximun wave heigh
print("maximun wave heigh: " + str(np.max(H)))
# mean period
print("mean period: " + str(np.mean(T)))
# significant wave period
print("significant wave period: " + str(highestN_stats(T, 3)))
# mean zero-crossing wave period
print("mean zero-crossing wave period: " + str(highestN_stats(T, 10)))
# maximun wave period
print("maximun wave period: " + str(np.max(T)))
mean wave heigh: 0.29362963917525775
root-mean-square wave height: 0.019087170539045063
significant wave heigh: 0.4431906976744186
mean of the highest one-tenth of waves: 0.5754684210526316
maximun wave heigh: 0.8243
mean period: 5.240359897172238
significant wave period: 7.038759689922486
mean zero-crossing wave period: 8.276315789473712
maximun wave period: 11.0
Wave Heigh Distribution¶
# Rayleigh distribution of individual wave heighs
fHs = [4.01 * (i/(4 * math.sqrt(moment(0, f, E)))**2) * np.exp(-2.005 * (i**2/(4 * np.sqrt(moment(0, f, E)))**2)) for i in H]
FHs = [1 - np.exp(-2.005 * (i**2/(4 * np.sqrt(moment(0, f, E)))**2)) for i in H]
# Plot - sort arrays
plt.figure(figsize=(5,5))
plt.hist(H, density = True, bins = 15, rwidth=0.95, color = "slateblue", alpha=0.5, label='Hist. Hi')
plt.plot(np.sort(H), np.array(fHs)[np.array(H).argsort()], c = "k", zorder = 10, label='Rayleight distribution')
plt.xlabel('H (m)')
plt.ylabel('f(H)')
plt.legend()
plt.xlim(np.nanmin(H), np.nanmax(H))
plt.show()
