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

\[S(f, \theta) = S(f)D(f, \theta).\]

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 and submodules. We start by importing numpy and MHKiT’s wave module.

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

buoy = '42012''swdir', buoy)
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 function. We see that the spectrum is calculated over ~1hour which results in 8572 different timestamps. We also see there are 47 frequency bins.

year = 2021
data_all =, year)
Dimensions:    (date: 8572, frequency: 47)
  * 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.1

Elevation 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 function. You can see the spectrum has the correct shape of 47 frequencies x 180 directions, and the units are indicated in the attributes.

date = np.datetime64('2021-02-21T12:40:00')
data = data_all.sel(date=date)
directions = np.arange(0, 360, 2.0)
spectrum =, directions)
<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]])
  * 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
    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 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.

<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]:, min=0.3)
<PolarAxesSubplot: title={'center': 'Elevation Variance Spectrum'}>
[7]:, min=0.3, fill=False, nlevels=4)
<PolarAxesSubplot: title={'center': 'Elevation Variance Spectrum'}>


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

[<matplotlib.lines.Line2D at 0x12e592670>]
spread =, directions), name="Spread", units="1")
<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.

rho = 1025 # kg/m^3
g = 9.81 # m/s^2*rho*g, name="Energy", units="J")
<PolarAxesSubplot: title={'center': 'Energy Spectrum'}>
[ ]: