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 dfarr: Frequency resolution DataArray.

property df

Frequency resolution numpy array.

TODO: Check if this can be removed in favor of dfarr.

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, 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).
rechunk (bool): Rechunk split dims so there is one single chunk.

Note:

Spectra are interpolated at fmin / fmax if they are not in self.freq.
Recommended rechunk==True so ufuncs with freq/dir as core dims will work.

to_energy(standard_name='sea_surface_wave_directional_energy_spectra')[source]: Convert from energy density (m2/Hz/degree) into wave energy spectra (m2).

hs(tail=True)[source]

Spectral significant wave height Hm0.

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 LH (2005). Waves in oceanic and coastal waters (page 82).

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.

momf(mom=0)[source]: Calculate given frequency moment.

momd(mom=0, theta=90.0)[source]: Calculate given directional moment.

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]

Directional wave spreading Dspr.

The one-sided directional width of the spectrum.

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). The statistical distribution of maxima of a random function. Proc. R. Soc. A237, 212-232.

sw(mask=nan)[source]

Spectral width parameter by Longuet-Higgins (1975).

Represents energy distribution over entire frequency range.

Args:

mask (float): value for missing data in output

Reference:

Longuet-Higgins (1975). On the joint distribution of the periods and amplitudes of sea waves. JGR, 80, 2688-2694.

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.

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.

partition(wsp_darr, wdir_darr, dep_darr, swells=3, agefac=1.7, wscut=0.3333)[source]

Partition wave spectra using Hanson’s watershed algorithm.

This method is not lazy, make sure array will fit into memory.

Args:

wsp_darr (DataArray): wind speed (m/s).
wdir_darr (DataArray): Wind direction (degree).
dep_darr (DataArray): Water depth (m).
swells (int): Number of swell partitions to compute.
agefac (float): Age factor.
wscut (float): Wind speed cutoff.

Returns:

part_spec (SpecArray): partitioned spectra with one extra dimension representig partition number.

Note:

Input DataArrays must have same non-spectral dims as SpecArray.

References:

Hanson, Jeffrey L., et al. “Pacific hindcast performance of three numerical wave models.” JTECH 26.8 (2009): 1614-1633.

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.

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=[0.01, 0.1, 1.0], 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”.

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)

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.

Note:

dataset needs to have lon/lat/time coordinates.
dataset with 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: - https://www.orcina.com/webhelp/OrcaFlex/Content/html/Wavetheory.htm - https://www.orcina.com/webhelp/OrcaFlex/Content/html/Directionconventions.htm

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)

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.

Note:

Only datasets with lat/lon coordinates are currently supported.
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.

Table of Contents

Wavespectra#

SpecArray#

SpecDataset#