# -*- coding: utf-8 -*-
"""
Functions for generating 1D and 2D wave spectra from parameters, modelled or observed data.
"""
#This module contains functions for calculating intergated wave parameters and
#reconstructing 2-D wave spectra from spectral parameters.
#
#Chris Old
#IES, School of Engineering, University of Edinburgh
#Aug 2021
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
# Non-Standard Python Dependencies
import numpy as np
from scipy.special import gamma as gammaFunc
# Local Module Dependencies
from .statsTools import circ_mod_pi
# Other Dependencies
# ----------------------------------------------------------------------------
# GLOBAL VARIABLES
# ----------------------------------------------------------------------------
grav = 9.81 # Mean gravitional acceleration at Earths surface m / s^2
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# FUNCTION DEFINITIONS
# ----------------------------------------------------------------------------
# ======== Theoretical 1-D Spectra Functions ============
[docs]
def pierson_moskowitz(f,Tp,Hs):
"""Generate Pierson-Moskowitz 1-D wave spectra
"""
try:
f = np.array(f)
except:
pass
assert isinstance(f, np.ndarray), 'f must be of type np.ndarray'
assert isinstance(Tp, (int,float)), 'Tp must be of type int or float'
assert isinstance(Hs, (int,float)), 'Hs must be of type int or float'
f.sort()
B_PM = (5/4)*(1/Tp)**4
A_PM = B_PM*(Hs/2)**2
Sf = A_PM*f**(-5)*np.exp(-B_PM*f**(-4))
return Sf
[docs]
def jonswap(f,Tp,Hs,gamma=None):
"""Generate Jonswap 1-D wave spectra
"""
try:
f = np.array(f)
except:
pass
assert isinstance(f, np.ndarray), 'f must be of type np.ndarray'
assert isinstance(Tp, (int,float)), 'Tp must be of type int or float'
assert isinstance(Hs, (int,float)), 'Hs must be of type int or float'
assert isinstance(gamma, (int,float, type(None))), \
'gamma must be of type int or float'
f.sort()
B_PM = (5/4)*(1/Tp)**4
A_PM = B_PM*(Hs/2)**2
S_f = A_PM*f**(-5)*np.exp(-B_PM*f**(-4))
if not gamma:
TpsqrtHs = Tp/np.sqrt(Hs);
if TpsqrtHs <= 3.6:
gamma = 5;
elif TpsqrtHs > 5:
gamma = 1;
else:
gamma = np.exp(5.75 - 1.15*TpsqrtHs);
# Cutoff frequencies for gamma function
siga = 0.07
sigb = 0.09
fp = 1/Tp # peak frequency
lind = np.where(f<=fp)
hind = np.where(f>fp)
Gf = np.zeros(f.shape)
Gf[lind] = gamma**np.exp(-(f[lind]-fp)**2/(2*siga**2*fp**2))
Gf[hind] = gamma**np.exp(-(f[hind]-fp)**2/(2*sigb**2*fp**2))
C = 1- 0.287*np.log(gamma)
Sf = C*S_f*Gf
return Sf
[docs]
def dsfunc1d(theta,thetap,dspr=None):
"""Generate 1-D wave directional spreading function
"""
ds1d = np.zeros(np.shape(theta))
dtheta = theta-thetap
dt = circ_mod_pi(np.deg2rad(dtheta))
if dspr is None:
indx = np.where(np.abs(dt) <= np.deg2rad(90.0))
ds1d[indx] = (2.0/np.pi)*np.square(np.cos(dt[indx]))
else:
s = (2.0/np.square(np.deg2rad(dspr)))-1.0
G1 = gammaFunc(s+1.0)
G2 = gammaFunc(s+0.5)
if G1 == np.Inf:
G1 = 1.0
G2 = 1.0
A2 = G1/(2.0*np.sqrt(np.pi)*G2)
ds1d[:] = A2*(np.cos(0.5*dt)**(2.0*s))[:]
return ds1d
# ============== 2-D Spectra Functions ==================
[docs]
def significant_wave_height(moments):
"""
Integrated Wave Parameter: Significant wave height
:param moments: Spectral moments calculated from the wave energy spectrum using :func:`WaveVal.WaveSpectra.spectralMoment()`.
:return: Hm0.
"""
Hm0 = 4.0 * np.sqrt(moments['m0'])
return Hm0
[docs]
def mean_wave_period(moments):
"""
Integrated Wave Parameter: Mean wave period
:param moments: Spectral moments calculated from the wave energy spectrum using :func:`WaveVal.WaveSpectra.spectralMoment()`.
:return: Tm01.
"""
Tm01 = (moments['m1']/moments['m0'])
return Tm01
[docs]
def zero_crossing_period(moments):
"""
Integrated Wave Parameter: Zero crossing period
:param moments: Spectral moments calculated from the wave energy spectrum using :func:`WaveVal.WaveSpectra.spectralMoment()`.
:return: Tm02.
"""
Tm02 = np.sqrt(moments['m0']/moments['m2'])
return Tm02
[docs]
def wave_energy_period(moments):
"""
Integrated Wave Parameter: Wave energy period
:param moments: Spectral moments calculated from the wave energy spectrum using :func:`WaveVal.WaveSpectra.spectralMoment()`.
:returns: Tm10.
"""
Tm10 = (moments['m0']/moments['m1'])
return Tm10
[docs]
def spectral_period(moments):
"""
Integrated Wave Parameter: Peak spectral period
:param moments: Spectral moments calculated from the wave energy spectrum using :func:`WaveVal.WaveSpectra.spectralMoment()`.
:returns: Tm02.
"""
Tdw2 = np.sqrt(moments['m-1']/moments['m1'])
return Tdw2
[docs]
def wave_steepness(moments):
"""
Integrated Wave Parameter: Wave steepness
:param moments: Spectral moments calculated from the wave energy spectrum using :func:`WaveVal.WaveSpectra.spectralMoment()`.
:returns: Tdw2.
"""
Hm0 = significant_wave_height(moments)
Tm01 = mean_wave_period(moments)
xi = (2.0 * np.pi / grav) * np.sqrt( Hm0 / np.square(Tm01) )
return xi
[docs]
def spectral_width(moments):
"""
Integrated Wave Parameter: Spectral width
:param moments: Spectral moments calculated from the wave energy spectrum using :func:`WaveVal.WaveSpectra.spectralMoment()`.
:returns: eps.
"""
num = moments['m0']*moments['m4'] - np.square(moments['m2'])
den = moments['m0']*moments['m4']
eps = np.sqrt(num / den)
return eps
[docs]
def spectral_narrowness(moments):
"""
Integrated Wave Parameter: Spectral narrowness
:param moments: Spectral moments calculated from the wave energy spectrum using :func:`WaveVal.WaveSpectra.spectralMoment()`.
:returns: nu.
"""
num = moments['m0']*moments['m2'] - np.square(moments['m1'])
den = np.square(moments['m1'])
nu = np.sqrt(num / den)
return nu
[docs]
def spectralMoment(f,df,S,moment):
"""
Calculate wave spectral moment
:param f: frequency bin centre values
:param df: frequency bin widths
:param S: spectral energy density in frquency bins
:param moment: interger representing moment to be calculated
:return: specMom
"""
if S.ndim > 1:
specMom = np.sum(pow(f,moment)*S*df,1)
else:
specMom = np.sum(pow(f,moment)*S*df)
return specMom
[docs]
def getSpectralMoments(f,df,S):
"""Get a set wave spectra moments.
"""
moments = {}
moms = [-1,0,1,2,4]
for mom in moms:
momStr = 'm'+str(mom)
moments[momStr] = spectralMoment(f,df,S,mom)
return moments
[docs]
def getSpectralParameters(moments):
"""Get wave spectra parameters for a set of moments.
"""
spec_params = {}
spec_params['Hm0'] = significant_wave_height(moments)
spec_params['Tm01'] = mean_wave_period(moments)
spec_params['Tm02'] = zero_crossing_period(moments)
spec_params['Tm10'] = wave_energy_period(moments)
spec_params['Tdw2'] = spectral_period(moments)
spec_params['wave_steepness'] = wave_steepness(moments)
spec_params['spectral_width'] = spectral_width(moments)
spec_params['spectral_narrowness'] = spectral_narrowness(moments)
return spec_params
# ============== 2-D Spectra Functions ==================
[docs]
def wavehgtvar(freq, df, dt, spec2D):
"""
Convert 2-D wave spectrum into wave height variance map in (f,theta) space.
"""
nfreq = spec2D.shape[0]
ntheta = spec2D.shape[1]
f_arry = np.outer(freq,np.ones(ntheta,))
df_arry = np.outer(df,np.ones(ntheta,))
dt_arry = np.outer(np.ones(nfreq,),np.radians(dt))
whv = spec2D*f_arry*df_arry*dt_arry
return whv
[docs]
def dsfuncnorm(freq, df, th1, sth1, ntheta):
"""
Normalised 2-D directional spreading function
Based on description given on NOAA `NDBC <https://www.ndbc.noaa.gov/measdes.shtml>`_ web site.
"""
nfreq = len(freq) # number of frequency bins
dth = 360.0 / ntheta # theta bin width
# Convert direction convention from meteorological to mathematical
theta1 = 270.0 - th1
dsFnc = np.empty((nfreq, ntheta), dtype=float)
# Apply weighting define in Earle et al, 1999 to avoid -ve energy
for ith in range(ntheta):
theta = ith * dth
dsFnc[:,ith] = np.exp(-0.5*np.square((np.radians(theta)-np.radians(theta1)+2.0*np.pi)/np.radians(sth1)))/(np.sqrt(2.0*np.pi)*np.radians(sth1))
dsFnc[:,ith] += np.exp(-0.5*np.square((np.radians(theta)-np.radians(theta1))/np.radians(sth1)))/(np.sqrt(2.0*np.pi)*np.radians(sth1))
dsFnc[:,ith] += np.exp(-0.5*np.square((np.radians(theta)-np.radians(theta1)-2.0*np.pi)/np.radians(sth1)))/(np.sqrt(2.0*np.pi)*np.radians(sth1))
# Normalise across theta
for ifrq in range(nfreq):
dsFnc[ifrq,:] = dsFnc[ifrq,:] / np.nansum(dsFnc[ifrq,:])
return dsFnc
[docs]
def spec2dnorm(freq, sf, df, th1, sth1, ntheta):
"""
Generate normalised 2-D directional spectrum from standard parameters retrieve
from a wavebuoy.
"""
# construct 2-D directional spreading function
dsFnc = dsfuncnorm(freq, df, th1, sth1, ntheta)
# construct 2-D directional spectrum
sp = np.outer(sf,np.ones(ntheta,))
sp2d = sp * dsFnc
return sp2d, dsFnc
[docs]
def dsfunc(freq, df, th1, th2, sth1, sth2, ntheta):
"""
Directional spreading function
Based on description of `MEDS
<https://forge.ifremer.fr/plugins/mediawiki/wiki/ww3/index.php/En:buoy_data:meds>`_
buoy data format given on the IFREMER Wiki pages.
:param freq:
:param df:
:param th1:
:param th2:
:param sth1:
:param sth2:
:param ntheta:
:return: dsFnc
"""
nfreq = len(freq) # number of frequency bins
dth = 360.0 / ntheta # theta bin width
W1 = 2.0/3.0
W2 = 1.0/6.0
if th2 is None:
th2 = np.zeros(th1.shape, dtype=float)
sth2 = np.zeros(sth1.shape, dtype = float)
W1 = 1.0
W2 = 0.0
# Convert direction convention from meteorological to mathematical
theta1 = 270.0 - th1
theta2 = 270.0 - th2
# Set theta2 to give minimum angular separation with theta1
th2mth1 = np.abs(theta2-theta1)
theta2[th2mth1 > 180.] = theta2[th2mth1 > 180.] + 180.
R1 = np.abs(1.0 - 0.5*(np.square(np.radians(sth1))))
R2 = np.abs(1.0 - 0.5*(np.square(np.radians(sth2))))
dsFnc = np.empty((nfreq, ntheta), dtype=float)
# Apply weighting define in Earle et al, 1999 to avoid -ve energy
for ith in range(ntheta):
theta = ith * dth
#dsFnc[:,ith] = ( 0.5 + R1*np.cos(np.radians(theta-theta1))
# + R2*np.cos(2.0*np.radians(theta-theta2)) ) / np.pi
#dsFnc[:,ith] = (dsFnc[:,ith]+np.abs(dsFnc[:,ith]))/2.0
dsFnc[:,ith] = ( 0.5 + W1*R1*np.cos(np.radians(theta-theta1))
+ W2*R2*np.cos(2.0*np.radians(theta-theta2)) ) / np.pi
# Normalise across theta
for ifrq in range(nfreq):
dsFnc[ifrq,:] = dsFnc[ifrq,:] / np.nansum(dsFnc[ifrq,:])
return dsFnc
[docs]
def spec2d(freq, sf, df, th1, th2, sth1, sth2, ntheta):
"""
Generate 2-D directional spectrum from standard parameters retrieve from
a wavebuoy.
"""
# construct 2-D directional spreading function
dsFnc = dsfunc(freq, df, th1, th2, sth1, sth2, ntheta)
# construct 2-D directional spectrum
sp = np.outer(sf,np.ones(ntheta,))
sp2d = sp * dsFnc
return sp2d, dsFnc
[docs]
def spec2dparam(freq,theta,thetap,Tp,Hs,gamma=None,dspr=None):
"""
Generate 2-D directional spectrum from parameters.
"""
ds = dsfunc1d(theta,thetap,dspr)
dsFnc = np.outer(np.ones(len(freq),),ds)
sf = jonswap(freq,Tp,Hs,gamma)
sp2d = np.outer(sf,np.ones(len(theta),))*dsFnc
return sp2d, dsFnc