Wavespectra#

Python library for ocean wave spectra.

Wavespectra is a library for processing ocean wave spectral data. It provides reading and writing of common spectral data formats, calculation of common integrated wave paramaters, several methods of spectral partitioning and spectral manipulation in a package focussed on speed and efficiency for large numbers of spectra.

The library built on top of xarray. The extended functionality is based on two accessor classes:

  • SpecArray: extends xarray’s DataArray with methods to manipulate wave spectra and calculate integrated spectral parameters.

  • SpecDataset: extends xarray’s Dataset with methods for writing spectra in different formats.

These accessors are registered to DataArray and Dataset objects as a new namespece called spec. To attach the namespace simply import the accessors into your code, e.g.:

import xarray as xr
import numpy as np
import datetime

from wavespectra.specarray import SpecArray

coords = {'time': [datetime.datetime(2017,01,n+1) for n in range(2)],
          'freq': [0.05,0.1],
          'dir': np.arange(0,360,120)}
efth = xr.DataArray(data=np.random.rand(2,2,3),
                    coords=coords,
                    dims=('time','freq', 'dir'),
                    name='efth')

In [1]: efth.spec.hs()
Out[1]:
<xarray.DataArray 'hs' (time: 2)>
array([ 10.128485,   9.510618])
Coordinates:
* time     (time) datetime64[ns] 2017-01-01 2017-01-02
Attributes:
    standard_name: sea_surface_wave_significant_height
    units: m

Using accessors in this way provides an elegant way to extend and leverage the power of xarray. In common usage, this mechanism is largely transparent, with the provided reading functions , returning the above classes from which the spectral methods can be directly accessed (see Input). For example, to read a swan spectral file and calculate the corresponing significant wave heights:

from wavespectra import read_swan

dset = read_swan('my_swan_file')
hs = dset.spec.hs()

Details of each class are below.

SpecArray#

SpecArray extends DataArray.

SpecArray provides several methods for calculating integrated wave parameters from wave spectra. For instance, significant wave height can be calculated from a hypothetical SpecArray object efth as:

hs = efth.spec.hs()

which returns a DataArray object hs.

The following attributes are required when defining SpecArray object:

  • Wave frequency coordinate named as freq, defined in \(Hz\).

  • Wave direction coordinate (if 2D spectra) named as dir, defined in \(degree\), (coming_from).

  • Wave energy density array named as efth, defined in:

    • 2D spectra \(E(\sigma,\theta)\): \(m^{2}{degree^{-1}}{s}\).

    • 1D spectra \(E(\sigma)\): \(m^{2}{s}\).

  • Time coordinate (if present) named as time.

The following methods are available in SpecArray:

class wavespectra.specarray.SpecArray(xarray_obj)[source]

Extends DataArray with methods to deal with wave spectra.

property freq

Frequency DataArray.

property dir

Direction DataArray.

property df

Frequency resolution DataArray.

property dd

Direction resolution float.

oned(skipna=True)[source]

Returns the one-dimensional frequency spectra.

Direction dimension is dropped after integrating.

Args:

  • skipna (bool): choose it to skip nans when integrating spectra. This is the default behaviour for sum() in DataArray. Notice it converts masks, where the entire array is nan, into zero.

split(fmin=None, fmax=None, dmin=None, dmax=None, interpolate=True, rechunk=True)[source]

Split spectra over freq and/or dir dims.

Args:

  • fmin (float): lowest frequency to split spectra, by default the lowest.

  • fmax (float): highest frequency to split spectra, by default the highest.

  • dmin (float): lowest direction to split spectra at, by default min(dir).

  • dmax (float): highest direction to split spectra at, by default max(dir).

  • interpolate (bool): Interpolate spectra at frequency cutoffs they are not available in coordinates of input spectra.

  • rechunk (bool): Rechunk split dims so there is one single chunk.

to_energy()[source]

Convert spectra from energy density (m2/Hz/degree) into wave energy (m2).

hs(tail=True)[source]

Spectral significant wave height Hm0.

Args:

  • tail (bool): if True fit high-frequency tail before integrating spectra.

hrms(tail=True)[source]

Root mean square wave height Hrms.

Args:

  • tail (bool): if True fit high-frequency tail before integrating spectra.

hmax()[source]

Maximum wave height Hmax.

hmax is the most probably value of the maximum individual wave height for each sea state. Note that maximum wave height can be higher (but not by much since the probability density function is rather narrow).

Reference:

  • Holthuijsen (2005).

scale_by_hs(expr, hs_min=-inf, hs_max=inf, tp_min=-inf, tp_max=inf, dpm_min=-inf, dpm_max=inf)[source]

Scale spectra using expression based on Significant Wave Height hs.

Args:

  • expr (str): expression to apply, e.g. ‘0.13*hs + 0.02’.

  • hs_min, hs_max, tp_min, tp_max, dpm_min, dpm_max (float): Ranges of hs,

    tp and dpm over which the scaling defined by expr is applied.

tp(smooth=True)[source]

Peak wave period Tp.

Args:

  • smooth (bool): True for the smooth wave period, False for the discrete period corresponding to the maxima in the frequency spectra.

Note:

  • The smooth wave period is defined from a parabolic fit around the discrete peak of the frequency spectrum.

fp(smooth=True)[source]

Peak wave frequency Fp.

Args:

  • smooth (bool): True for the smooth wave period, False for the discrete period corresponding to the maxima in the frequency spectrum.

Note:

  • The smooth wave period is defined from a parabolic fit around the discrete peak of the frequency spectrum.

momf(mom=0)[source]

Frequency moment.

Args:

  • mom (int): Moment to calculate.

Returns:

  • mf (DataArray): The mth frequency moments for each direction.

momd(mom=0, theta=90.0)[source]

Directional moment.

Args:

  • mom (int): Moment to calculate.

  • theta (float): angle offset.

Returns:

  • msin (DataArray): Sin component of the mth directional moment for each frequency.

  • mcos (DataArray): Cosine component of the mth directional moment for each frequency.

tm01()[source]

Mean absolute wave period Tm01.

True average period from the 1st spectral moment.

tm02()[source]

Mean absolute wave period Tm02.

Average period of zero up-crossings (Zhang, 2011).

dm()[source]

Mean wave direction from the 1st spectral moment Dm.

dp()[source]

Peak wave direction Dp.

Defined as the direction where the energy density of the frequency-integrated spectrum is maximum.

dpm()[source]

Peak wave direction Dpm.

Note From WW3 Manual:

  • peak wave direction, defined like the mean direction, using the frequency/wavenumber bin containing of the spectrum F(k) that contains the peak frequency only.

dspr()[source]

Mean directional wave spread Dspr.

The one-sided directional width of the spectrum.

fdspr(mom=1)[source]

Directional wave spread at frequency \(dspr(f)\).

The directional width of the spectrum at each frequency.

Args:

  • mom (int): Directional moment for calculating the mth directional spread.

dpspr(mom=1)[source]

Peak directional wave spread Dpspr.

The directional width of the spectrum at peak frequency.

Args:

  • mom (int): Directional moment to calculate the mth directional spread.

crsd(theta=90.0)[source]

Add description.

swe()[source]

Spectral width parameter by Cartwright and Longuet-Higgins (1956).

Represents the range of frequencies where the dominant energy exists.

Reference:

  • Cartwright and Longuet-Higgins (1956).

sw()[source]

Spectral width parameter by Longuet-Higgins (1975).

Represents energy distribution over entire frequency range.

Reference:

  • Longuet-Higgins (1975).

gw()[source]

Gaussian frequency spread by Bunney et al. (2014).

Represents gaussian width of a swell partition.

Reference:

  • Bunney et al. (2014).

alpha(smooth=True)[source]

Jonswap fetch dependant scaling coefficient alpha.

Args:

  • smooth (bool): True for the smooth wave frequency, False for the discrete frequency corresponding to the peak in the frequency spectrum.

Reference:

  • Phillips (1957).

gamma(smooth=True, scaled=True)[source]

Jonswap peak enhancement factor gamma.

Represents the ratio between the peak in the frequency spectrum \(E(f)\) and its associate Pierson-Moskowitz shape.

Args:

  • smooth (bool): True for the smooth wave frequency, False for the discrete frequency corresponding to the peak in the frequency spectrum.

  • scaled (bool): True to apply polynomial approximation to gamma.

goda()[source]

Goda peakedness parameter.

Reference:

  • Goda (1970).

celerity(depth=None)[source]

Wave celerity C from frequency coords.

Args:

  • depth (float): Water depth, use deep water approximation by default.

Returns;

  • C: ndarray of same shape as freq with wave celerity for each frequency.

uss_x(depth=None, theta=90.0)[source]

Stokes drift - x component, at sea surface. No high frequency tail.

Args:

  • depth (float): Water depth, use deep water approximation by default.

  • theta (float): angle offset.

uss_y(depth=None, theta=90.0)[source]

Stokes drift - y component, at sea surface. No high frequency tail.

Args:

  • depth (float): Water depth, use deep water approximation by default.

  • theta (float): angle offset.

uss(depth=None)[source]

Stokes drift - speed, at sea surface. No high frequency tail.

Args:

  • depth (float): Water depth, use deep water approximation by default.

mss(depth=None)[source]

Mean squared slope of sea surface.

Args:

  • depth (float): Water depth, use deep water approximation by default.

wavelen(depth=None)[source]

Wavelength L from frequency coords.

Args:

  • depth (float): Water depth, use deep water approximation by default.

Returns;

  • L: ndarray of same shape as freq with wavelength for each frequency.

stats(stats, fmin=None, fmax=None, dmin=None, dmax=None, names=None)[source]

Calculate multiple spectral stats into a Dataset.

Args:

  • stats (list): strings specifying stats to be calculated. (dict): keys are stats names, vals are dicts with kwargs to use with corresponding method.

  • fmin (float): lower frequencies for splitting spectra before calculating stats.

  • fmax (float): upper frequencies for splitting spectra before calculating stats.

  • dmin (float): lower directions for splitting spectra before calculating stats.

  • dmax (float): upper directions for splitting spectra before calculating stats.

  • names (list): strings to rename each stat in output Dataset.

Returns:

  • Dataset with all spectral statistics specified.

Note:

  • All stats names must correspond to methods implemented in this class.

  • If names is provided, its length must correspond to the length of stats.

rotate(angle) → DataArray[source]

Rotate spectra.

Args:

  • angle (float): Angle in degrees to rotate spectra, positive for clockwise and negative for counter-clockwise rotation.

Returns:

  • efth (DataArray): Rotated spectra.

Note:

  • Directions are interpolated so that the rotated spectra have the same directional bins as the original spectra.

smooth(freq_window=3, dir_window=3)[source]

Smooth spectra with a running average.

Args:

  • freq_window (int): Rolling window size along freq dim.

  • dir_window (int): Rolling window size along dir dim.

Returns:

  • efth (DataArray): Smoothed spectra.

Note:

  • Window sizes must be odd to ensure symmetry.

interp(freq=None, dir=None, maintain_m0=True)[source]

Interpolate onto new spectral basis.

Args:

  • freq (DataArray, 1darray): Frequencies of interpolated spectra (Hz).

  • dir (DataArray, 1darray): Directions of interpolated spectra (deg).

  • maintain_m0 (bool): Ensure variance is conserved in interpolated spectra.

Returns:

  • dsi (DataArray): Regridded spectra.

Note:

  • All freq below lowest freq are interpolated assuming \(E_d(f=0)=0\).

  • \(Ed(f)\) is set to zero for new freq above the highest freq in dset.

  • Only the ‘linear’ method is currently supported.

interp_like(other, maintain_m0=True)[source]

Interpolate onto coordinates from other spectra.

Args:

  • other (Dataset, DataArray): Spectra defining new spectral basis.

  • maintain_m0 (bool): Ensure variance is conserved in interpolated spectra.

Returns:

  • dsi (DataArray): Regridded spectra.

plot(kind='contourf', normalised=True, logradius=True, as_period=False, rmin=None, rmax=None, show_theta_labels=True, show_radii_labels=True, radii_ticks=None, radii_labels_angle=22.5, radii_labels_size=8, cbar_ticks=None, cmap='RdBu_r', extend='neither', efth_min=0.001, **kwargs)[source]

Plot spectra in polar axis.

Args:

  • kind (str): Plot kind, one of (contourf, contour, pcolormesh).

  • normalised (bool): Show efth normalised between 0 and 1.

  • logradius (bool): Set log radii.

  • as_period (bool): Set radii as wave period instead of frequency.

  • rmin (float): Minimum value to clip the radius axis.

  • rmax (float): Maximum value to clip the radius axis.

  • show_theta_labels (bool): Show direction tick labels.

  • show_radii_labels (bool): Show radii tick labels.

  • radii_ticks (array): Tick values for radii.

  • radii_labels_angle (float): Polar angle at which radii labels are positioned.

  • radii_labels_size (float): Fontsize for radii labels.

  • cbar_ticks (array): Tick values for colorbar.

  • cmap (str, obj): Colormap to use.

  • efth_min (float): Clip energy density below this value.

  • kwargs: All extra kwargs are passed to the plotting method defined by kind.

Returns:

  • pobj: The xarray object returned by calling da.plot.{kind}(**kwargs).

Note:

  • If normalised==True, contourf uses a logarithmic colour scale by default.

  • Plot and axes can be redefined from the returned xarray object.

  • Xarray uses the sharex, sharey args to control which panels receive axis labels. In order to set labels for all panels, set these to False.

  • Masking of low values can be done in contourf by setting efth_min larger than the lowest contour level along with extend set to “neither” or “min”.

rmse(other)[source]

Root-Mean-Square Error among spectral bins.

Args:

  • other (SpecArray): Spectra to calculate binwise RMSE against.

Note:

  • Spectral coordinates in self and other must be the same.

  • Non-spectral coordinates in self and other should ideally be the same. Coordinates are attempted to be broadcast if they are different.

Reference:

  • Bunney et al. (2014).

fit_jonswap(gamma0=1.5, spectra=True, params=True)[source]

Nonlinear fit Jonswap spectra.

Args:

  • gamma0 (float): Initial guess for gamma.

  • spectra (bool): Return fitted spectra.

  • params (bool): Return fitted parameters.

Return:

  • dsout (Dataset): Fitted Jonswap spectra and/or Jonswap parameters peak frequency, significant wave height and peak enhancement factor gamma.

Note:

  • At least one of spectra or params must be True.

fit_gaussian(gw0=1.5, spectra=True, params=True)[source]

Nonlinear fit Gaussian width.

Args:

  • gw0 (float): Initial guess for gw.

  • spectra (bool): Return fitted spectra.

  • params (bool): Return fitted parameters.

Return:

  • dsout (Dataset): Fitted Gaussian spectra and/or Gaussian parameters peak frequency, significant wave height and gaussian width.

Note:

  • At least one of spectra or params must be True.

SpecDataset#

SpecDataset extends Dataset.

SpecDataset is useful for storing wave spectra alongside other variables that share some common dimensions. It provides methods for writing wave spectra into different file formats.

SpecDataset works as a wrapper around SpecArray. All public methods from SpecArray can be directly accessed from SpecDataset. For instance, these two calls are equivalent:

# Calling hs method from SpecArray
hs = dset.efth.spec.hs()
# Calling hs method from SpecDataset
hs = dset.spec.hs()

both cases return identical DataArray objects hs.

The following methods are available in SpecDataset:

class wavespectra.specdataset.SpecDataset(xarray_dset)[source]

Extends xarray’s Dataset to deal with wave spectra datasets.

Plugin functions defined in wavespectra/output/<module> are attached as methods in this accessor class.

sel(lons, lats, method='idw', tolerance=2.0, dset_lons=None, dset_lats=None, **kwargs)[source]

Select stations near or at locations defined by (lons, lats) vector.

Args:

  • lons (list): Longitude values of locations to select.

  • lats (list): Latitude values of locations to select.

  • method (str): Method to use for inexact matches:

    • idw: Inverse distance weighting selection.

    • nearest: Nearest site selection.

    • bbox: Sites inside bbox [min(lons), min(lats)], [max(lons), max(lats)].

    • None: Only exact matches.

  • tolerance (float): Maximum distance between locations and original stations for inexact matches.

  • dset_lons (array): Longitude of stations in dset, not required but could help improove speed.

  • dset_lats (array): Latitude of stations in dset, not required but could help improove speed.

  • kwargs: Extra keywargs to pass to the respective sel function (i.e., sel_nearest, sel_idw).

Return:

  • dset (SpecDataset): Stations Dataset selected at locations defined by zip(lons, lats).

Note:

  • tolerance behaves differently with methods ‘idw’ and ‘nearest’. In ‘idw’ sites with no neighbours within tolerance are masked whereas in ‘nearest’ an exception is raised.

  • dset_lons, dset_lats are not required but can improve performance when dset is chunked with site=1 (expensive to access site coords) and improve precision if projected coors are provided at high latitudes.

to_funwave(filename, clip=True)

Write spectra in Funwave format.

Args:

  • filename (str): Name for output Funwave file.

  • clip (bool): Clip directions outside [-90, 90] range in cartesian convention.

Note:

  • Format description: https://fengyanshi.github.io/build/html/wavemaker_para.html.

  • Directions converted from wavespectra (0N, CW, from) to Cartesian (0E, CCW, to).

  • Funwave only seems to deal with directions in the [-90, 90] range in cartesian convention, use clip=True to clip spectra outside that range.

  • Both 2D \(E(f,d)\) and 1d \(E(f)\) spectra are supported.

  • If the SpecArray has more than one spectrum, multiple files are created in a zip archive defined by replacing the extension of filename by “.zip”.

to_json(filename, mode='w', date_format='%Y-%m-%dT%H:%M:%SZ')

Write spectra in json format.

Xarray’s to_dict it used to dump dataset into dictionary to save as a json file.

Args:

  • filename (str): name of output json file.

  • mode (str): file mode, by default w (create or overwrite).

  • date_format(str): strftime format for serializing datetimes.

to_netcdf(filename, specname='efth', ncformat='NETCDF4', compress=True, packed=True, time_encoding={'units': 'seconds since 1970-01-01'})

Write spectra in netCDF format using wavespectra conventions.

Args:

  • filename (str): name of output netcdf file.

  • specname (str): name of spectra variable in dataset.

  • ncformat (str): netcdf format for output, see options in native to_netcdf method.

  • compress (bool): if True output is compressed, has no effect for NETCDF3.

  • packed (bool): Pack spectra as int32 dtype and 1e-5 scale_factor.

  • time_encoding (dict): force standard time units in output files.

to_octopus(filename, site_id='spec', fcut=0.125, missing_val=-99999, ntime=None, lons=None, lats=None, compresslevel=6)

Save spectra in Octopus format.

Args:

  • filename (str): name for output OCTOPUS file.

  • site_id (str): used to construct LPoint header.

  • fcut (float): frequency for splitting spectra.

  • missing_value (int): missing value in output file.

  • ntime (int, None): number of times to load into memory before dumping output file if full dataset does not fit into memory, choose None to load all times.

  • lons: (np.array, None): longitudes to use for each site if not in dataset.

  • lats: (np.array, None): latitudes to use for each site if not in dataset.

  • compresslevel (int): compression level for gzip compression (1-9).

Note:

  • Output files are gzipped if filename ends with .gz.

  • lons/lats parameters may be prescribed to set site locations if lon/lat are not variables in the dataset (their sizes must match the number of sites).

  • If lons/lats are not specified and the dataset does not have lon/lat coords, all coordinates default to zero.

  • multiple locations dumped at same file with one location header per site.

  • 1D spectra not supported.

  • ntime=None optimises speed as the dataset is loaded into memory however the dataset may not fit into memory in which case a smaller number of times may be prescribed.

to_orcaflex(model, minEnergy=1e-06)

Writes the spectrum to an Orcaflex model

Uses the orcaflex API (OrcFxAPI) to set the wave-data of the provided orcaflex model.

The axis system conversion used is:

  • Orcaflex global X = Towards East

  • Orcaflex global Y = Towards North

This function creates a wave-train for each of the directions in this spectrum using a user-defined spectrum.

Calculation of wave-components in orcaflex is computationally expensive. To save computational time:

  1. Use the minEnergy parameter of this function to define a treshold for the amount of energy in a wave-train.

  2. In orcaflex itself: limit the amount of wave-components.

  3. Before exporting: regrid the spectrum to a lower amount of directions.

Orcaflex theory:

Example:
>>> from OrcFxAPI import *
>>> from wavespectra import read_triaxys
>>> m = Model()
>>> spectrum = read_triaxys("triaxys.DIRSPEC")).isel(time=0)  # get only the fist spectrum in time
>>> spectrum.spec.to_orcaflex(m)
Args:

  • model : orcaflex model (OrcFxAPI.model instance).

  • minEnergy [1e-6] : threshold for minimum sum of energy in a direction before it is exported.

Note:

  • an Orcaflex license is required to work with the orcaflex API.

  • Only 2D spectra \(E(f,d)\) are currently supported.

  • The DataArray should contain only a single spectrum. Hint: first_spetrum = spectra.isel(time=0).

to_swan(filename, append=False, id='Created by wavespectra', ntime=None, lons=None, lats=None, compresslevel=6)

Write spectra in SWAN ASCII format.

Args:

  • filename (str): str, name for output SWAN ASCII file.

  • append (bool): if True append to existing filename.

  • id (str): used for header in output file.

  • ntime (int, None): number of times to load into memory before dumping output file if full dataset does not fit into memory, choose None to load all times.

  • lons: (np.array, None): longitudes to use for each site, if None use.

  • lats: (np.array, None): latitudes to use for each site, if None use.

  • compresslevel (int): compression level for gzip compression (1-9).

Note:

  • Output files are gzipped if filename ends with .gz.

  • lons/lats parameters may be prescribed to set site locations if lon/lat are not variables in the dataset (their sizes must match the number of sites).

  • If lons/lats are not specified and the dataset does not have lon/lat coords, all coordinates default to zero.

  • Extra dimensions other than time, site, lon, lat, freq, dim not yet supported.

  • Only 2D spectra E(f,d) are currently supported.

  • ntime=None optimises speed as the dataset is loaded into memory however the dataset may not fit into memory in which case a smaller number of times may be prescribed.

to_ww3(filename, ncformat='NETCDF4', compress=False)

Save spectra in native WW3 netCDF format.

Args:

  • filename (str): name of output WW3 netcdf file.

  • ncformat (str): netcdf format for output, see options in native to_netcdf method.

  • compress (bool): if True output is compressed, has no effect for NETCDF3.

On this page