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.
Reading spectra from files#
Several methods are provided to read various file formats including spectral wave models like WW3, SWAN, WWM and observation instruments such as TRIAXYS and SPOTTER.
In [1]: import matplotlib.pyplot as plt
In [2]: from wavespectra import read_ww3
In [3]: dset = read_ww3("_static/ww3file.nc")
In [4]: dset
Out[4]:
<xarray.Dataset>
Dimensions: (time: 9, site: 2, freq: 25, dir: 24)
Coordinates:
* freq (freq) float32 0.04118 0.0453 0.04983 ... 0.3352 0.3687 0.4056
* site (site) int32 1 2
* time (time) datetime64[ns] 2014-12-01 2014-12-01T12:00:00 ... 2014-12-05
* dir (dir) float32 270.0 255.0 240.0 225.0 ... 330.0 315.0 300.0 285.0
Data variables:
dpt (time, site) float32 dask.array<chunksize=(9, 2), meta=np.ndarray>
efth (time, site, freq, dir) float32 dask.array<chunksize=(9, 2, 25, 24), meta=np.ndarray>
lat (site) float32 dask.array<chunksize=(2,), meta=np.ndarray>
lon (site) float32 dask.array<chunksize=(2,), meta=np.ndarray>
wspd (time, site) float32 dask.array<chunksize=(9, 2), meta=np.ndarray>
wdir (time, site) float32 dask.array<chunksize=(9, 2), meta=np.ndarray>
The spec namespace#
Wavespectra defines a new namespace accessor called spec which is attached to xarray objects. This namespace provides access to methods from the two main objects in wavespectra:
which extend functionality from xarray’s DataArray and Dataset objects respectively.
SpecArray#
In [5]: dset.efth.spec
Out[5]:
<SpecArray 'efth' (time: 9, site: 2, freq: 25, dir: 24)>
dask.array<SpecArray shape=(9, 2, 25, 24), dtype=float32, chunksize=(9, 2, 25, 24), chunktype=numpy.ndarray>
Coordinates:
* freq (freq) float32 0.04118 0.0453 0.04983 ... 0.3352 0.3687 0.4056
* site (site) int32 1 2
* time (time) datetime64[ns] 2014-12-01 2014-12-01T12:00:00 ... 2014-12-05
* dir (dir) float32 270.0 255.0 240.0 225.0 ... 330.0 315.0 300.0 285.0
Attributes:
standard_name: sea_surface_wave_directional_variance_spectral_density
units: m2 s degree-1
SpecDset#
In [6]: dset.spec
Out[6]:
<SpecDataset>
Dimensions: (time: 9, site: 2, freq: 25, dir: 24)
Coordinates:
* freq (freq) float32 0.04118 0.0453 0.04983 ... 0.3352 0.3687 0.4056
* site (site) int32 1 2
* time (time) datetime64[ns] 2014-12-01 2014-12-01T12:00:00 ... 2014-12-05
* dir (dir) float32 270.0 255.0 240.0 225.0 ... 330.0 315.0 300.0 285.0
Data variables:
dpt (time, site) float32 dask.array<SpecDataset>
efth (time, site, freq, dir) float32 dask.array<SpecDataset>
lat (site) float32 dask.array<SpecDataset>
lon (site) float32 dask.array<SpecDataset>
wspd (time, site) float32 dask.array<SpecDataset>
wdir (time, site) float32 dask.array<SpecDataset>
Spectral wave parameters#
Several methods are available to calculate integrated wave parameters. They can be accessed from both SpecArray (efth variable) and SpecDset accessors:
In [7]: hs = dset.efth.spec.hs()
In [8]: hs
Out[8]:
<xarray.DataArray 'hs' (time: 9, site: 2)>
dask.array<mul, shape=(9, 2), dtype=float64, chunksize=(9, 2), chunktype=numpy.ndarray>
Coordinates:
* site (site) int32 1 2
* time (time) datetime64[ns] 2014-12-01 2014-12-01T12:00:00 ... 2014-12-05
Attributes:
standard_name: sea_surface_wave_significant_height
units: m
In [9]: hs1 = dset.spec.hs()
In [10]: hs.identical(hs1)
Out[10]: True
In [11]: hs.plot.line(x="time");
In [12]: plt.draw()
In [13]: stats = dset.spec.stats(
....: ["hs", "hmax", "tp", "tm01", "tm02", "dpm", "dm", "dspr", "swe"]
....: )
....:
In [14]: stats
Out[14]:
<xarray.Dataset>
Dimensions: (site: 2, time: 9)
Coordinates:
* site (site) int32 1 2
* time (time) datetime64[ns] 2014-12-01 2014-12-01T12:00:00 ... 2014-12-05
Data variables:
hs (time, site) float64 dask.array<chunksize=(9, 2), meta=np.ndarray>
hmax (time, site) float64 dask.array<chunksize=(9, 2), meta=np.ndarray>
tp (time, site) float32 dask.array<chunksize=(9, 2), meta=np.ndarray>
tm01 (time, site) float64 dask.array<chunksize=(9, 2), meta=np.ndarray>
tm02 (time, site) float64 dask.array<chunksize=(9, 2), meta=np.ndarray>
dpm (time, site) float32 dask.array<chunksize=(9, 2), meta=np.ndarray>
dm (time, site) float32 dask.array<chunksize=(9, 2), meta=np.ndarray>
dspr (time, site) float64 dask.array<chunksize=(9, 2), meta=np.ndarray>
swe (time, site) float64 dask.array<chunksize=(9, 2), meta=np.ndarray>
Attributes:
standard_name: sea_surface_wave_significant_height
units: m
In [15]: fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(8, 6))
In [16]: stats.hs.plot.line(ax=ax1, x="time");
In [17]: stats.hmax.plot.line(ax=ax2, x="time");
In [18]: stats.dpm.plot.line(ax=ax3, x="time");
In [19]: stats.dspr.plot.line(ax=ax4, x="time");
In [20]: stats.tp.plot.line(ax=ax5, x="time");
In [21]: stats.tm01.plot.line(ax=ax6, x="time");
In [22]: plt.draw()
Partitioning#
Two different partitioning techniques are available, a simple spectral split based on frequency / direction cutoffs and the watershed algorithm of Hanson et al. (2008).
Spectral split#
In [23]: fcut = 1 / 8
In [24]: sea = dset.spec.split(fmin=fcut)
In [25]: swell = dset.spec.split(fmax=fcut)
In [26]: dset.freq.values
Out[26]:
array([0.04118 , 0.045298 , 0.0498278 , 0.05481058, 0.06029164,
0.06632081, 0.07295289, 0.08024818, 0.08827299, 0.09710029,
0.10681032, 0.11749136, 0.1292405 , 0.14216454, 0.15638101,
0.17201911, 0.18922101, 0.20814312, 0.22895744, 0.25185317,
0.27703848, 0.30474234, 0.3352166 , 0.36873826, 0.40561208],
dtype=float32)
In [27]: sea.freq.values
Out[27]:
array([0.125 , 0.1292405 , 0.14216454, 0.15638101, 0.17201911,
0.18922101, 0.20814312, 0.22895744, 0.25185317, 0.27703848,
0.30474234, 0.33521661, 0.36873826, 0.40561208])
In [28]: swell.freq.values
Out[28]:
array([0.04118 , 0.045298 , 0.0498278 , 0.05481058, 0.06029164,
0.06632081, 0.07295289, 0.08024818, 0.08827299, 0.09710029,
0.10681032, 0.11749136, 0.125 ])
In [29]: dset.spec.hs().isel(site=0).plot(label='Full spectrum', marker='o');
In [30]: sea.spec.hs().isel(site=0).plot(label='Sea', marker='o');
In [31]: swell.spec.hs().isel(site=0).plot(label='Swell', marker='o');
In [32]: plt.draw()
Watershed#
In [33]: dspart = dset.spec.partition(dset.wspd, dset.wdir, dset.dpt)
In [34]: pstats = dspart.spec.stats(["hs", "dpm"])
In [35]: pstats
Out[35]:
<xarray.Dataset>
Dimensions: (site: 2, time: 9, part: 4)
Coordinates:
* site (site) int32 1 2
* time (time) datetime64[ns] 2014-12-01 2014-12-01T12:00:00 ... 2014-12-05
* part (part) int64 0 1 2 3
Data variables:
hs (part, time, site) float64 dask.array<chunksize=(4, 9, 2), meta=np.ndarray>
dpm (part, time, site) float32 dask.array<chunksize=(4, 9, 2), meta=np.ndarray>
Attributes:
standard_name: sea_surface_wave_significant_height
units: m
In [36]: fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8))
In [37]: hs.isel(site=0).plot(ax=ax1, label='Full spectrum', marker='o');
In [38]: pstats.hs.isel(part=0, site=0).plot(ax=ax1, label='Partition 0 (sea)', marker='o');
In [39]: pstats.hs.isel(part=1, site=0).plot(ax=ax1, label='Partition 1 (swell 1)', marker='o');
In [40]: pstats.hs.isel(part=2, site=0).plot(ax=ax1, label='Partition 2 (swell 2)', marker='o');
In [41]: pstats.hs.isel(part=3, site=0).plot(ax=ax1, label='Partition 3 (swell 3)', marker='o');
In [42]: dset.spec.dpm().isel(site=0).plot(ax=ax2, label='Full spectrum', marker='o');
In [43]: pstats.dpm.isel(part=0, site=0).plot(ax=ax2, label='Partition 0 (sea)', marker='o');
In [44]: pstats.dpm.isel(part=1, site=0).plot(ax=ax2, label='Partition 1 (swell 1)', marker='o');
In [45]: pstats.dpm.isel(part=2, site=0).plot(ax=ax2, label='Partition 2 (swell 2)', marker='o');
In [46]: pstats.dpm.isel(part=3, site=0).plot(ax=ax2, label='Partition 3 (swell 3)', marker='o');
In [47]: plt.draw()
Plotting#
Wavespectra wraps the plotting functionality from xarray to allow easily defining frequency-direction spectral plots in polar coordinates.
In [48]: ds = dset.isel(site=0, time=[0, 1]).spec.split(fmin=0.05, fmax=2.0)
In [49]: 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 [50]: import cmocean
In [51]: ds = dset.isel(site=0).spec.split(fmax=0.18).spec.oned().rename({"freq": "period"})
In [52]: ds = ds.assign_coords({"period": 1 / ds.period})
In [53]: ds.period.attrs.update({"standard_name": "sea_surface_wave_period", "units": "s"})
In [54]: ds.plot.contourf(x="time", y="period", vmax=1.25, cmap=cmocean.cm.thermal, levels=10);
Selecting#
Wavespectra complements xarray’s selecting and interpolating functionality with functions to select and
interpolate from site coordinates with the sel method.
In [55]: idw = dset.spec.sel(
....: lons=[92.01, 92.05, 92.09],
....: lats=[19.812, 19.875, 19.935],
....: method="idw"
....: )
....:
In [56]: idw
Out[56]:
<xarray.Dataset>
Dimensions: (freq: 25, time: 9, dir: 24, site: 3)
Coordinates:
* freq (freq) float32 0.04118 0.0453 0.04983 ... 0.3352 0.3687 0.4056
* time (time) datetime64[ns] 2014-12-01 2014-12-01T12:00:00 ... 2014-12-05
* dir (dir) float32 270.0 255.0 240.0 225.0 ... 330.0 315.0 300.0 285.0
* site (site) int64 0 1 2
Data variables:
dpt (time, site) float32 dask.array<chunksize=(9, 1), meta=np.ndarray>
efth (time, site, freq, dir) float32 dask.array<chunksize=(9, 1, 25, 24), meta=np.ndarray>
lat (site) float64 19.81 19.88 19.93
lon (site) float64 92.01 92.05 92.09
wspd (time, site) float32 dask.array<chunksize=(9, 1), meta=np.ndarray>
wdir (time, site) float32 dask.array<chunksize=(9, 1), meta=np.ndarray>
In [57]: p = plt.scatter(dset.lon, dset.lat, 80, dset.isel(time=0).spec.hs(), label="Dataset points");
In [58]: p = plt.scatter(idw.lon, idw.lat, 100, idw.isel(time=0).spec.hs(), marker="v", label="Interpolated point");
In [59]: plt.draw()
The nearest neighbour and bbox options are also available besides inverse distance weighting (idw).