Directional Wave Spectrum from NDBC Buoys
Some NDBC buoys measure the full directional wave spectrum. A directional wave spectrum \(S(f,\theta)\) contains the sea elevation information over all different directions. This can be used, for instance, to (1) reconstruct a sea surface elevation time-history over some area (as opposed to at just one point), or (2) to calculate the effect of wave loads on a non-axissymetric structure. It is constructed from the omnidirectional spectrum \(S(f)\) and a spread function \(D(f, \theta)\) as
NDBC stores this information in 5 parameters, each stored in a different file. One parameter (swden) is simply the omnidirectional spectrum, while the other 4 (swdir, swdir2, swr1, swr2) are used to reconstruct the spread function. For more information on the the parameters and reconstruction see the Spectral Wave Data section in this NDBC site.
In this example we will: * check if a specific buoy has directional information and request that data * reconstruct and plot the directional wave spectrum \(S(f, \theta)\) * take a deeper look at the omnidirectional spectrum \(S(f)\) and spread function \(D(f, \theta)\) separately * calculate and plot the energy spectrum
We use functions from both the wave.io.ndbc and wave.graphics submodules. We start by importing numpy and MHKiT’s wave module.
[1]:
import numpy as np
from mhkit import wave
Request Data
First, we will check if the buoy we’re interested in–42012—has directional data. To do this we use the wave.io.ndbc.available_data function with any one of the 4 parameters used for the spread function. We can see that this buoy does have directional data for the years 2009-2021.
[2]:
buoy = '42012'
wave.io.ndbc.available_data('swdir', buoy)
[2]:
| id | year | filename | |
|---|---|---|---|
| 312 | 42012 | 2009 | 42012d2009.txt.gz |
| 313 | 42012 | 2010 | 42012d2010.txt.gz |
| 314 | 42012 | 2011 | 42012d2011.txt.gz |
| 315 | 42012 | 2012 | 42012d2012.txt.gz |
| 316 | 42012 | 2013 | 42012d2013.txt.gz |
| 317 | 42012 | 2014 | 42012d2014.txt.gz |
| 318 | 42012 | 2015 | 42012d2015.txt.gz |
| 319 | 42012 | 2016 | 42012d2016.txt.gz |
| 320 | 42012 | 2017 | 42012d2017.txt.gz |
| 321 | 42012 | 2018 | 42012d2018.txt.gz |
| 322 | 42012 | 2019 | 42012d2019.txt.gz |
| 323 | 42012 | 2020 | 42012d2020.txt.gz |
| 324 | 42012 | 2021 | 42012d2021.txt.gz |
We will ask for the data for 2021. To do this we use the wave.io.ndbc.request_directional_data function. We see that the spectrum is calculated over ~1hour which results in 8572 different timestamps. We also see there are 47 frequency bins.
[3]:
year = 2021
data_all = wave.io.ndbc.request_directional_data(buoy, year)
data_all
[3]:
<xarray.Dataset>
Dimensions: (date: 8572, frequency: 47)
Coordinates:
* date (date) datetime64[ns] 2021-01-01T00:40:00 ... 2021-12-31T23:40:00
* frequency (frequency) float64 0.02 0.0325 0.0375 ... 0.445 0.465 0.485
Data variables:
swden (date, frequency) float64 0.0 0.0 0.0 0.0 ... 0.01 0.01 0.02 0.0
swdir (date, frequency) float64 194.0 13.0 6.0 ... 173.0 167.0 159.0
swdir2 (date, frequency) float64 191.0 11.0 9.0 ... 178.0 169.0 106.0
swr1 (date, frequency) float64 0.16 0.26 0.21 0.41 ... 0.76 0.87 0.63
swr2 (date, frequency) float64 0.99 0.93 0.9 0.84 ... 0.28 0.67 0.1Elevation Variance Spectrum
We now use the data we requested and reconstruct the directional elevation variance spectrum. The first step is to decide on the resolution of wave direction. Here we calculate the spectrum every 2 degrees. For the reconstruction we use the wave.io.ndbc.create_directional_spectrum function. You can see the spectrum has the correct shape of 47 frequencies x 180 directions, and the units are indicated in the attributes.
[4]:
date = np.datetime64('2021-02-21T12:40:00')
data = data_all.sel(date=date)
directions = np.arange(0, 360, 2.0)
spectrum = wave.io.ndbc.create_directional_spectrum(data, directions)
spectrum
[4]:
<xarray.DataArray (frequency: 47, direction: 180)>
array([[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
...,
[2.21512935e-04, 3.00007246e-04, 3.89537948e-04, ...,
4.78972678e-05, 9.59294370e-05, 1.53638415e-04],
[5.05171904e-04, 5.68137407e-04, 6.34119449e-04, ...,
3.34637188e-04, 3.88407564e-04, 4.45255232e-04],
[2.26590428e-04, 2.84552135e-04, 3.45964024e-04, ...,
7.35307830e-05, 1.21080223e-04, 1.72097824e-04]])
Coordinates:
* frequency (frequency) float64 0.02 0.0325 0.0375 ... 0.445 0.465 0.485
* direction (direction) float64 0.0 2.0 4.0 6.0 ... 352.0 354.0 356.0 358.0
Attributes:
units: m^2/Hz/deg
long_name: Elevation variance spectrum
standard_name: spectrum
description: *Elevation variance (m^2)* spectrum (/Hz/deg).We finally plot the spectrum using the wave.graphics.plot_directional_spectrum function. Notice that there are negative values, which is due to NDBC’s analysis approach which results in negative values far from the mean direction. From this NDBC site: > \(D(f,\theta)\) can take on negative values because of the trigonometric sine and cosine functions. There are several approaches to prevent or deal with the negative values. For more
information and discussion of some approaches see: Use of advanced directional wave spectra analysis methods, M. D. Earle, K. E. Steele, and D. W. C. Wang, Ocean Engineering, Volume 26, Issue 12, December 1999, Pages 1421-1434.
[5]:
wave.graphics.plot_directional_spectrum(spectrum)
[5]:
<PolarAxesSubplot: title={'center': 'Elevation Variance Spectrum'}>
We now show a few other ways of plotting this directional spectra. In the first plot below we set a minimum value of \(0.3m^2/Hz\). In the second plot, we plot 3 contour lines with no fill.
[6]:
wave.graphics.plot_directional_spectrum(spectrum, min=0.3)
[6]:
<PolarAxesSubplot: title={'center': 'Elevation Variance Spectrum'}>
[7]:
wave.graphics.plot_directional_spectrum(spectrum, min=0.3, fill=False, nlevels=4)
[7]:
<PolarAxesSubplot: title={'center': 'Elevation Variance Spectrum'}>
Components
In this section we will look plot both the omnidirectional spectrum \(S(f)\) and the spread function \(D(f, \theta)\) separately. Remember that the spectrum above was given by \(S(f)D(f, \theta)\).
[8]:
data['swden'].plot()
[8]:
[<matplotlib.lines.Line2D at 0x12e592670>]
[9]:
spread = wave.io.ndbc.create_spread_function(data, directions)
wave.graphics.plot_directional_spectrum(spread, name="Spread", units="1")
[9]:
<PolarAxesSubplot: title={'center': 'Spread Spectrum'}>
Energy Spectrum
Finally, we show how to calculate and plot other spectra. Here we show the energy spectrum, which is proportional to the elevation variance spectrum with a proportionality constant of \(\rho g\). That is, \(S_E(f, \theta) = \rho g S(f, \theta)\). We provide the plotting function additional information so that it displays the correct name and units.
[10]:
rho = 1025 # kg/m^3
g = 9.81 # m/s^2
wave.graphics.plot_directional_spectrum(spectrum*rho*g, name="Energy", units="J")
[10]:
<PolarAxesSubplot: title={'center': 'Energy Spectrum'}>
[ ]: