Dfsu - Speed and direction

import numpy as np
import mikeio

ds = mikeio.read("../tests/testdata/HD2D.dfsu")
ds

<mikeio.Dataset>
Geometry: Dfsu2D
Dimensions: (time:9, element:884)
Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00 (9 records)
Items:
  0:  Surface elevation <Surface Elevation> (meter)
  1:  U velocity <u velocity component> (meter per sec)
  2:  V velocity <v velocity component> (meter per sec)
  3:  Current speed <Current Speed> (meter per sec)

This file is missing current direction :-(

Lets’fix that!

Calculate speed & direction

ds.U_velocity

<mikeio.DataArray>
Name: U velocity
Geometry: Dfsu2D
Dimensions: (time:9, element:884)
Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00 (9 records)

In order to use Numpy functions on a DataArray, we first convert the DataArrays (U, V) to standard NumPy ndarrays.

u = ds.U_velocity.to_numpy()
v = ds.V_velocity.to_numpy()

direction = np.mod(90 -np.rad2deg(np.arctan2(v,u)),360)

Write new file

from mikeio.eum import ItemInfo, EUMUnit, EUMType

ds["Current direction"] = mikeio.DataArray(direction, time= ds.time, item = ItemInfo("Current direction", EUMType.Current_Direction, EUMUnit.degree), geometry=ds.geometry)
ds

<mikeio.Dataset>
Geometry: Dfsu2D
Dimensions: (time:9, element:884)
Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00 (9 records)
Items:
  0:  Surface elevation <Surface Elevation> (meter)
  1:  U velocity <u velocity component> (meter per sec)
  2:  V velocity <v velocity component> (meter per sec)
  3:  Current speed <Current Speed> (meter per sec)
  4:  Current direction <Current Direction> (degree)

ds.to_dfs("speed_direction.dfsu")

nds = mikeio.read("speed_direction.dfsu")
nds

<mikeio.Dataset>
Geometry: Dfsu2D
Dimensions: (time:9, element:884)
Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00 (9 records)
Items:
  0:  Surface elevation <Surface Elevation> (meter)
  1:  U velocity <u velocity component> (meter per sec)
  2:  V velocity <v velocity component> (meter per sec)
  3:  Current speed <Current Speed> (meter per sec)
  4:  Current direction <Current Direction> (degree)

Plot

step = 1

ax = ds.Current_speed[step].plot(figsize=(10,10))
ax.set_ylim([None, 6903000])
ax.set_xlim([607000, None])
ec = ds.geometry.element_coordinates
x = ec[:,0]
y = ec[:,1]
u = ds.U_velocity.to_numpy()
v = ds.V_velocity.to_numpy()
ax.quiver(x, y, u[step], v[step], scale=6, minshaft=3);

Plot quiver on Cartesian overlay instead

Create overset grid and interpolate data on to this

g = ds.geometry.get_overset_grid(dx=50)
g

<mikeio.Grid2D>
x-axis: nx=33 points from x0=605902 to x1=607502 with dx=50
y-axis: ny=93 points from y0=6.90242e+06 to y1=6.90702e+06 with dy=50
Number of grid points: 3069

g.projection

'UTM-29'

ui = ds.U_velocity.interp_like(g)
vi = ds.V_velocity.interp_like(g)

ui.plot();

ax = ds.Current_speed.plot(figsize=(10,10))

u = ui.to_numpy()
v = vi.to_numpy()
ax.quiver(g.x, g.y, u[step], v[step], scale=8, minshaft=5)

ax.set_ylim([None, 6903000])
ax.set_xlim([607000, None])
ax.set_title(f'Current speed with overset grid; {ds.time[step]}')
ax.set_xlabel("Easting (m)")
ax.set_ylabel("Northing (m)");

Clean up

import os
os.remove("speed_direction.dfsu")