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: (dir: 24, time: 9, site: 2, freq: 25)
Coordinates:
* dir (dir) float32 270.0 255.0 240.0 225.0 ... 330.0 315.0 300.0 285.0
* 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
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:
* dir (dir) float32 270.0 255.0 240.0 225.0 ... 330.0 315.0 300.0 285.0
* 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
Attributes:
standard_name: sea_surface_wave_directional_variance_spectral_density
units: m2 s degree-1
SpecDset
In [6]: dset.spec
Out[6]:
<SpecDataset>
Dimensions: (dir: 24, time: 9, site: 2, freq: 25)
Coordinates:
* dir (dir) float32 270.0 255.0 240.0 225.0 ... 330.0 315.0 300.0 285.0
* 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
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: (dir: 24, freq: 25, time: 9, site: 3)
Coordinates:
* dir (dir) float32 270.0 255.0 240.0 225.0 ... 330.0 315.0 300.0 285.0
* 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
* 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).