MIKE 21 SW can output full spectral information in points, along lines or in an area. In all these cases data are stored in dfsu files with additional axes: frequency and directions.
This notebook explores reading full spectral dfsu files from MIKE 21 SW as
import numpy as np
import matplotlib.pyplot as plt
from mikeio import Dfsu
fn = "../tests/testdata/pt_spectra.dfsu"
dfs = Dfsu(fn)
dfs
DfsuSpectral0D Number of frequencies: 25 Number of directions: 16 Projection: LONG/LAT Items: 0: Point 1: Energy density <Wave energy density> (meter pow 2 sec per deg) Time: 31 steps with dt=600.0s 2017-10-27 00:00:00 -- 2017-10-27 05:00:00
ds = dfs.read(time_steps=0)
spec = np.squeeze(ds[0])
ax = dfs.plot_spectrum(spec, rmax=8, plot_type="patch");
dird = np.round(dfs.directions*180/np.pi,2)
ax.set_thetagrids(dird,labels=dird);
Data in dfsu line spectra is node-based contrary to must other dfsu-formats.
fn = "../tests/testdata/line_spectra.dfsu"
dfs = Dfsu(fn)
dfs
DfsuSpectral1D Number of nodes: 10 Number of frequencies: 25 Number of directions: 16 Projection: LONG/LAT Items: 0: Energy density <Wave energy density> (meter pow 2 sec per deg) Time: 4 steps with dt=6000.0s 2017-10-27 00:00:00 -- 2017-10-27 05:00:00
ds = dfs.read()
spec = np.squeeze(ds[0][0,3,:,:]) # note first 3 points are outside domain
dfs.plot_spectrum(spec, cmap="Greys", rmax=8, r_as_periods=True);
Hm0 = dfs.calc_Hm0_from_spectrum(ds[0])
timestep = 0
plt.plot(dfs.node_coordinates[:,0],Hm0[timestep,:])
plt.title('Hm0 on a line crossing the domain')
plt.xlabel("Longitude [degrees]")
plt.ylabel("Hm0 [m]");
fn = "../tests/testdata/area_spectra.dfsu"
dfs = Dfsu(fn)
dfs
DfsuSpectral2D Number of elements: 40 Number of nodes: 33 Number of frequencies: 25 Number of directions: 16 Projection: LONG/LAT Items: 0: Energy density <Wave energy density> (meter pow 2 sec per deg) Time: 3 steps with dt=9000.0s 2017-10-27 00:00:00 -- 2017-10-27 05:00:00
ds = dfs.read()
Hm0 = dfs.calc_Hm0_from_spectrum(ds[0])
Hm0.shape
(3, 40)
ax = dfs.plot(Hm0[-1,:], label="Hm0 [m]") # last time step
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude");
from ipywidgets import interact
from datetime import timedelta
@interact
def plot_element(id=(0,dfs.n_elements-1), step=(0,dfs.n_timesteps-1)):
spec = np.squeeze(ds[0][step,id])
time = dfs.start_time + timedelta(seconds=(step*dfs.timestep))
dfs.plot_spectrum(spec, vmax=0.04, vmin=0, rmax=8, title=f"Wave spectrum, {time}, element: {id}");
plt.show()
interactive(children=(IntSlider(value=19, description='id', max=39), IntSlider(value=1, description='step', ma…