Quick start#
Wavespectra is an open source project for processing ocean wave spectral data. The library is built on top of xarray and provides reading and writing of different spectral data formats, calculation of common integrated wave paramaters, spectral partitioning and spectral manipulation in a package focussed on speed and efficiency for large numbers of spectra.
Wavespectra works by defining some conventions for spectral data in xarray objects:
In [1]: import numpy as np
In [2]: import xarray as xr
In [3]: import wavespectra
# Spectral data in m2/Hz/degree
In [4]: efth = np.array([0.0000e+00, 1.0942e-03, 5.8909e-02, 2.0724e-01,
...: 2.0267e-01, 1.6886e-01, 1.1980e-01, 7.7867e-02,
...: 4.9664e-02, 3.2361e-02, 2.0902e-02, 1.3622e-02,
...: 8.8645e-03, 5.6254e-03, 3.3722e-03, 2.0276e-03,
...: 1.1350e-03, 4.5993e-04, 2.0441e-04, 0.0000e+00])
...:
# Frequency in Hz
In [5]: freq = np.array([0.05 , 0.075, 0.1 , 0.125, 0.15 , 0.175, 0.2 ,
...: 0.225, 0.25 , 0.275, 0.3 , 0.325, 0.35 , 0.375,
...: 0.4 , 0.425, 0.45 , 0.475, 0.5 , 0.525])
...:
# Direction in degree
In [6]: dir = np.array([0])
# The DataArray and Dataset objects below are ready to use with wavespectra
In [7]: da = xr.DataArray(
...: data=np.expand_dims(efth, 1),
...: dims=["freq", "dir"],
...: coords=dict(freq=freq, dir=dir),
...: name="efth",
...: )
...:
In [8]: da.spec.hs()
Out[8]:
<xarray.DataArray 'hs' ()> Size: 8B
array(0.62439675)
Attributes:
standard_name: sea_surface_wave_significant_height
units: m
In [9]: dset = da.to_dataset()
In [10]: dset.spec.hs()
Out[10]:
<xarray.DataArray 'hs' ()> Size: 8B
array(0.62439675)
Attributes:
standard_name: sea_surface_wave_significant_height
units: m
The spec namespace#
Once wavespectra has been imported, DataArray and Dataset objects following this
convention will have the spec accessor available, from which the functionality of
wavespectra can be accessed.
SpecArray#
The SpecArray accessor extends DataArray objects
and is the main entry point for spectral data manipulation methods.
In [11]: dset.efth.spec
Out[11]:
<SpecArray 'efth' (freq: 20, dir: 1)> Size: 160B
array([[0.0000e+00],
[1.0942e-03],
[5.8909e-02],
[2.0724e-01],
[2.0267e-01],
[1.6886e-01],
[1.1980e-01],
[7.7867e-02],
[4.9664e-02],
[3.2361e-02],
[2.0902e-02],
[1.3622e-02],
[8.8645e-03],
[5.6254e-03],
[3.3722e-03],
[2.0276e-03],
[1.1350e-03],
[4.5993e-04],
[2.0441e-04],
[0.0000e+00]])
Coordinates:
* freq (freq) float64 160B 0.05 0.075 0.1 0.125 ... 0.45 0.475 0.5 0.525
* dir (dir) int64 8B 0
SpecDataset#
The SpecDataset accessor extends Dataset objects
to allow accessing SpecArray methods directly from
the Dataset, in addition to providing other methods such as writing spectral data files.
In [12]: dset.spec
Out[12]:
<SpecDataset> Size: 328B
Dimensions: (freq: 20, dir: 1)
Coordinates:
* freq (freq) float64 160B 0.05 0.075 0.1 0.125 ... 0.45 0.475 0.5 0.525
* dir (dir) int64 8B 0
Data variables:
efth (freq, dir) float64 160B 0.0 0.001094 0.05891 ... 0.0002044 0.0
Reading spectra from files#
Wavespectra provides functions to read various file formats including spectral wave models like WAVEWATCHIII, SWAN and WWM, observation instruments such as TRIAXYS and SPOTTER, and industry standard formats including ERA5, NDBC, Octopus among others.
In [13]: import matplotlib.pyplot as plt
In [14]: from wavespectra import read_ww3
In [15]: dset = read_ww3("_static/ww3file.nc")
In [16]: dset
Out[16]:
<xarray.Dataset> Size: 44kB
Dimensions: (time: 9, site: 2, freq: 25, dir: 24)
Coordinates:
* freq (freq) float32 100B 0.04118 0.0453 0.04983 ... 0.3352 0.3687 0.4056
* site (site) int32 8B 1 2
* time (time) datetime64[ns] 72B 2014-12-01 ... 2014-12-05
* dir (dir) float32 96B 270.0 255.0 240.0 225.0 ... 315.0 300.0 285.0
Data variables:
dpt (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
efth (time, site, freq, dir) float32 43kB dask.array<chunksize=(9, 2, 25, 24), meta=np.ndarray>
lat (site) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
lon (site) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
wspd (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
wdir (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
In version 4, xarray engines have been defined for all wavespectra readers, allowing for direct reading of spectral data using xarray.open_dataset.
In [17]: dset = xr.open_dataset("_static/ww3file.nc", engine="ww3")
Spectral wave parameters#
Several methods are available to calculate integrated wave parameters. They can be accessed from both SpecArray (efth variable) and SpecDataset accessors:
In [18]: hs = dset.efth.spec.hs()
In [19]: hs
Out[19]:
<xarray.DataArray 'hs' (time: 9, site: 2)> Size: 72B
dask.array<mul, shape=(9, 2), dtype=float32, chunksize=(9, 2), chunktype=numpy.ndarray>
Coordinates:
* site (site) int32 8B 1 2
* time (time) datetime64[ns] 72B 2014-12-01 ... 2014-12-05
Attributes:
standard_name: sea_surface_wave_significant_height
units: m
In [20]: hs1 = dset.spec.hs()
In [21]: hs.identical(hs1)
Out[21]: True
In [22]: hs.plot.line(x="time");
In [23]: plt.draw()
In [24]: stats = dset.spec.stats(
....: ["hs", "hmax", "tp", "tm01", "tm02", "dpm", "dm", "dspr", "swe"]
....: )
....:
In [25]: stats
Out[25]:
<xarray.Dataset> Size: 800B
Dimensions: (site: 2, time: 9)
Coordinates:
* site (site) int32 8B 1 2
* time (time) datetime64[ns] 72B 2014-12-01 ... 2014-12-05
Data variables:
hs (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
hmax (time, site) float64 144B dask.array<chunksize=(9, 2), meta=np.ndarray>
tp (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
tm01 (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
tm02 (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
dpm (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
dm (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
dspr (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
swe (time, site) float32 72B dask.array<chunksize=(9, 2), meta=np.ndarray>
Attributes:
standard_name: sea_surface_wave_significant_height
units: m
In [26]: fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(8, 6))
In [27]: stats.hs.plot.line(ax=ax1, x="time");
In [28]: stats.hmax.plot.line(ax=ax2, x="time");
In [29]: stats.dpm.plot.line(ax=ax3, x="time");
In [30]: stats.dspr.plot.line(ax=ax4, x="time");
In [31]: stats.tp.plot.line(ax=ax5, x="time");
In [32]: stats.tm01.plot.line(ax=ax6, x="time");
In [33]: plt.draw()
Interpolate#
A custom interpolation method takes care of the cyclic nature of the wave direction.
In [34]: ds = dset.efth.isel(site=0, time=0).sortby("dir")
In [35]: freq = np.arange(ds.freq.min(), ds.freq.max()+0.001, 0.001)
In [36]: dir = np.arange(0, 360, 1)
In [37]: ds_interp = ds.spec.interp(freq=freq, dir=dir)
In [38]: fig, axs = plt.subplots(1, 2, figsize=(12, 5))
In [39]: ds.plot(ax=axs[0], x="dir", y="freq", cmap="turbo", add_colorbar=False)
Out[39]: <matplotlib.collections.QuadMesh at 0x7fea67c83130>
In [40]: ds_interp.plot(ax=axs[1], x="dir", y="freq", cmap="turbo", add_colorbar=False)
Out[40]: <matplotlib.collections.QuadMesh at 0x7fea67d54e50>
In [41]: plt.draw()
Smoothing#
Spectra smoothing is available using a running average method.
In [42]: ds_smooth = ds.spec.smooth(freq_window=5, dir_window=5)
In [43]: dss = xr.concat([np.log10(ds), np.log10(ds_smooth)], dim="smooth")
In [44]: dss["smooth"] = ["false", "true"]
In [45]: dss.plot(col="smooth", x="dir", y="freq", cmap="turbo", add_colorbar=False);
In [46]: plt.draw()
Spectra file writing#
Several methods are available in the SpecDataset accessor for writing spectral data to different file formats. The following example writes the dataset to a SWAN ASCII file:
In [47]: dset.spec.to_swan("specfile.swn")
In [48]: !head -n 40 specfile.swn
SWAN 1 Swan standard spectral file
$ Created by wavespectra
$
TIME time-dependent data
1 time coding option
LONLAT locations in spherical coordinates
2 number of locations
92.099998 19.950001
92.000000 19.799999
AFREQ absolute frequencies in Hz
25 number of frequencies
0.04118
0.04530
0.04983
0.05481
0.06029
0.06632
0.07295
0.08025
0.08827
0.09710
0.10681
0.11749
0.12924
0.14216
0.15638
0.17202
0.18922
0.20814
0.22896
0.25185
0.27704
0.30474
0.33522
0.36874
0.40561
NDIR spectral nautical directions in degr
24 number of directions
270.0000
255.0000
Plotting#
Wavespectra wraps the plotting functionality from xarray to allow easy plotting of frequency-direction spectral plots in polar coordinates.
In [49]: ds = dset.isel(site=0, time=[0, 1]).spec.split(fmin=0.05, fmax=0.4)
In [50]: ds.spec.plot(
....: kind="contourf",
....: col="time",
....: as_period=False,
....: normalised=True,
....: logradius=True,
....: add_colorbar=False,
....: figsize=(8, 5)
....: );
....:
Plotting Hovmoller diagrams of frequency spectra timeseries can be done in only a few lines.
In [51]: import cmocean
In [52]: ds = dset.isel(site=0).spec.split(fmax=0.18).spec.oned().rename({"freq": "period"})
In [53]: ds = ds.assign_coords({"period": 1 / ds.period})
In [54]: ds.period.attrs.update({"standard_name": "sea_surface_wave_period", "units": "s"})
In [55]: ds.plot.contourf(x="time", y="period", vmax=1.25, cmap=cmocean.cm.thermal, levels=10);
Partitioning#
Different partitioning techniques are available within the spec.partition namespace. The partitioning methods follow the naming convention defined in the WAVEWATCHIII model (ptm1, ptm2, etc) with the addition of some custom methods. In the following example, the ptm1 method is used to partition the dataset into wind sea and three swells (ptm1 is equivalent to the former spec.partition() method deprecated in version 4).
In [56]: dspart = dset.spec.partition.ptm1(dset.wspd, dset.wdir, dset.dpt)
In [57]: pstats = dspart.spec.stats(["hs", "dpm"])
In [58]: pstats
Out[58]:
<xarray.Dataset> Size: 688B
Dimensions: (site: 2, time: 9, part: 4)
Coordinates:
* site (site) int32 8B 1 2
* time (time) datetime64[ns] 72B 2014-12-01 ... 2014-12-05
* part (part) int64 32B 0 1 2 3
Data variables:
hs (part, time, site) float32 288B dask.array<chunksize=(4, 9, 2), meta=np.ndarray>
dpm (part, time, site) float32 288B dask.array<chunksize=(4, 9, 2), meta=np.ndarray>
Attributes:
standard_name: sea_surface_wave_significant_height
units: m
In [59]: fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8))
In [60]: hs.isel(site=0).plot(ax=ax1, label='Full spectrum', marker='o');
In [61]: pstats.hs.isel(part=0, site=0).plot(ax=ax1, label='Partition 0 (sea)', marker='o');
In [62]: pstats.hs.isel(part=1, site=0).plot(ax=ax1, label='Partition 1 (swell 1)', marker='o');
In [63]: pstats.hs.isel(part=2, site=0).plot(ax=ax1, label='Partition 2 (swell 2)', marker='o');
In [64]: pstats.hs.isel(part=3, site=0).plot(ax=ax1, label='Partition 3 (swell 3)', marker='o');
In [65]: dset.spec.dpm().isel(site=0).plot(ax=ax2, label='Full spectrum', marker='o');
In [66]: pstats.dpm.isel(part=0, site=0).plot(ax=ax2, label='Partition 0 (sea)', marker='o');
In [67]: pstats.dpm.isel(part=1, site=0).plot(ax=ax2, label='Partition 1 (swell 1)', marker='o');
In [68]: pstats.dpm.isel(part=2, site=0).plot(ax=ax2, label='Partition 2 (swell 2)', marker='o');
In [69]: pstats.dpm.isel(part=3, site=0).plot(ax=ax2, label='Partition 3 (swell 3)', marker='o');
In [70]: plt.draw()
Construction#
Spectral construction functionality has been implemented in version 4 with different shape functions available for frequency and direction such as Jonswap and Cartwright:
In [71]: from wavespectra.construct import construct_partition
In [72]: freq = np.arange(0.03, 0.401, 0.001)
In [73]: dir = np.arange(0, 360, 1)
In [74]: ds = construct_partition(
....: freq_name="jonswap",
....: dir_name="cartwright",
....: freq_kwargs={"freq": freq, "fp": 0.1, "gamma": 3.3, "hs": 1.5},
....: dir_kwargs={"dir": dir, "dm": 60, "dspr": 30},
....: )
....:
In [75]: ds.spec.plot();
In [76]: ds.spec.oned().plot(figsize=(8, 4));
Selecting#
Wavespectra complements xarray’s selecting and interpolating functionality with functions to select and
interpolate from site coordinates with the sel method.
In [77]: idw = dset.spec.sel(
....: lons=[92, 92.05, 92.1, 92.1, 92.1, 92.1, 92.05, 92, 92, 92],
....: lats=[19.8, 19.8, 19.8, 19.85, 19.9, 19.95, 19.95, 19.95, 19.9, 19.85],
....: method="idw"
....: )
....:
In [78]: idw
Out[78]:
<xarray.Dataset> Size: 218kB
Dimensions: (time: 9, site: 10, freq: 25, dir: 24)
Coordinates:
* freq (freq) float32 100B 0.04118 0.0453 0.04983 ... 0.3352 0.3687 0.4056
* time (time) datetime64[ns] 72B 2014-12-01 ... 2014-12-05
* dir (dir) float32 96B 270.0 255.0 240.0 225.0 ... 315.0 300.0 285.0
* site (site) int64 80B 0 1 2 3 4 5 6 7 8 9
Data variables:
dpt (time, site) float32 360B dask.array<chunksize=(9, 1), meta=np.ndarray>
efth (time, site, freq, dir) float32 216kB dask.array<chunksize=(9, 1, 25, 24), meta=np.ndarray>
lat (site) float64 80B 19.8 19.8 19.8 19.85 ... 19.95 19.95 19.9 19.85
lon (site) float64 80B 92.0 92.05 92.1 92.1 ... 92.05 92.0 92.0 92.0
wspd (time, site) float32 360B dask.array<chunksize=(9, 1), meta=np.ndarray>
wdir (time, site) float32 360B dask.array<chunksize=(9, 1), meta=np.ndarray>
In [79]: p = plt.scatter(dset.lon, dset.lat, s=250, c=dset.isel(time=0).spec.hs(), cmap="turbo", marker="$\u25EF$", label="Dataset points");
In [80]: p = plt.scatter(idw.lon, idw.lat, s=80, c=idw.isel(time=0).spec.hs(), cmap="turbo", marker="o", edgecolor="k", label="Interpolated point");
In [81]: plt.draw()
The nearest neighbour and bbox options are also available besides inverse distance weighting (idw).