Dfsu

Dfsu and mesh files are both flexible mesh file formats used by MIKE 21/3 engines. The .mesh file is an ASCII file for storing the flexible mesh geometry. The .dfsu file is a binary dfs file with data on this mesh. The mesh geometry is available in a .dfsu file as static items.

For a detailed description of the .mesh and .dfsu file specification see the flexible file format documentation.

The flexible mesh

The mesh geometry in a .mesh or a .dfsu file consists of a list of nodes and a list of elements.

Each node has:

  • Node id
  • x,y,z coordinates
  • Code (0 for internal water points, 1 for land, >1 for open boundary)

Each element has:

  • Element id
  • Element table; specifies for each element the nodes that defines the element. (the number of nodes defines the type: triangular, quadrilateral, prism etc.)
Note

In MIKE Zero, node ids, element ids and layer ids are 1-based. In MIKE IO, all ids are 0-based following standard Python indexing. That means, as an example, that when finding the element closest to a point its id will be 1 lower in MIKE IO compared to examining the file in MIKE Zero.

MIKE IO Flexible Mesh Geometry

MIKE IO has Flexible Mesh Geometry classes, e.g. GeometryFM2D, containing the list of node coordinates and the element table which defines the mesh, as well as a number of derived properties (e.g. element coordinates) and methods making it convenient to work with the mesh.

If a .dfsu file is read with mikeio.read, the returned Dataset ds will contain a Flexible Mesh Geometry geometry. If a .dfsu or a .mesh file is opened with mikeio.open, the returned object will also contain a Flexible Mesh Geometry geometry.

import mikeio

ds = mikeio.read("../data/oresundHD_run1.dfsu")
ds.geometry
Flexible Mesh Geometry: Dfsu2D
number of nodes: 2046
number of elements: 3612
projection: UTM-33
dfs = mikeio.open("../data/oresundHD_run1.dfsu")
dfs.geometry
Flexible Mesh Geometry: Dfsu2D
number of nodes: 2046
number of elements: 3612
projection: UTM-33

Common Dfsu and Mesh properties

MIKE IO has Dfsu classes for .dfsu files and a Mesh class for .mesh files which both have a mikeio.spatial.GeometryFM2D/mikeio.spatial.GeometryFM3D accessible through the ´geometry´ accessor.

Dfsu types

The following dfsu file types are supported by MIKE IO.

  • 2D horizontal.
  • 3D layered.
  • 2D vertical profile - a vertical slice through a 3D layered file.
  • 1D vertical column - a vertical dfs1 file and is produced by taking out one column of a 3D layered file.
  • 3D/4D SW, two horizontal dimensions and 1-2 spectral dimensions. Output from MIKE 21 SW.

When a dfsu file is opened with mikeio.open() the returned dfs object will be a specialized class Dfsu2DH, Dfsu3D, Dfsu2DV, or DfsuSpectral according to the type of dfsu file.

The layered files (3d, 2d/1d vertical) can have both sigma- and z-layers or only sigma-layers.

In most cases values are stored in cell centers and vertical (z) information in nodes, but the following values types exists:

  • Standard value type, storing values on elements and/or nodes. This is the default type.
  • Face value type, storing values on element faces. This is used e.g. for HD decoupling files, to store the discharge between elements.
  • Spectral value type, for each node or element, storing vales for a number of frequencies and/or directions. This is the file type for spectral output from the MIKE 21 SW.

Visualization

Geometry

The geometry (bathymetry) can be visualized in various ways.

geom  = dfs.geometry
geom.plot();

geom.plot.mesh();

geom.plot(vmin=-20);

geom.plot.contour(levels=[-20,-10], cmap='tab10');

geom.plot.contourf(levels=6, vmin=-30);

Data

The data can be visualized in a similar way.

da = ds["Surface elevation"]
da.plot(cmap='plasma');

ax = da.plot.contour(show_mesh=True, cmap="tab20")
ax.set_ylim(6135000, 6160000);

Subplots

Pass a matplotlib axes to create side-by-side comparisons:

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 7))

ds["Surface elevation"].isel(time=-1).plot.contourf(ax=ax1, title="Surface elevation")
ds["V velocity"].isel(time=-1).plot.contourf(
    ax=ax2, title="V velocity", label="Northward velocity [m/s]"
);

Working with data

Subsetting by element properties

Elements can be selected based on their coordinates, e.g. by depth:

import numpy as np

ze = dfs.geometry.element_coordinates[:, 2]
deep_elems = np.where(ze < -20)[0]

ds_deep = dfs.read(items=0, time=-1, elements=deep_elems)
ax = ds_deep[0].plot()
dfs.geometry.plot.mesh(ax=ax);

Boolean masking

Instead of reading a subset, you can mask values directly on a DataArray:

da_last = ds["Surface elevation"].isel(time=-1).copy()
da_last[ze >= -20] = np.nan
da_last.plot(show_mesh=True);

3D layered files: z-coordinates

Layered 3D dfsu files (and the vertical-profile / vertical-column slices derived from them) carry a time-varying vertical coordinate at every mesh node. MIKE IO exposes this through the .z accessor on a DataArray (or Dataset):

ds3d = mikeio.read("../data/oresund_sigma_z.dfsu")
da3d = ds3d["Temperature"]
da3d.z.nodes.shape       # (n_time, n_nodes) — z at every mesh node, every timestep
(3, 12042)
da3d.z.elements.shape    # (n_time, n_elements) — mean of each element's node z's
(3, 17118)

da.z.nodes returns the raw node z-coordinates exactly as stored in the file; da.z.elements is derived as the per-element mean of the column-node z’s and is cached on first access. The same accessor is available on Datasets (ds.z.nodes, ds.z.elements) and is delegated to the first item.

Worked example: dz-weighted column-mean

Element-center z-coordinates are useful for vertical integration. The cell height dz between layer interfaces can be approximated from successive element-center z-values within a column; for a clean dz, take a vertical column slice first:

column = ds3d.sel(x=333934.1, y=6158101.5)
ze = column["Temperature"].z.elements   # (n_time, n_layers)
dz = np.diff(ze, axis=-1)
dz = np.concatenate([dz[:, :1], dz], axis=-1)   # pad to match n_layers

T = column["Temperature"].to_numpy()
T_mean = (T * dz).sum(axis=-1) / dz.sum(axis=-1)
T_mean.shape   # depth-weighted temperature per timestep
(3,)

Non-layered DataArrays

For any DataArray whose geometry is not layered (dfs2, 2D dfsu, spectral, …), accessing .z returns a placeholder that raises AttributeError to make the geometry constraint discoverable:

ds2 = mikeio.read("../data/random.dfs2")
ds2[0].z.nodes
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[16], line 2
      1 ds2 = mikeio.read("../data/random.dfs2")
----> 2 ds2[0].z.nodes

File ~/work/mikeio/mikeio/src/mikeio/dataset/_z_accessor.py:119, in NullZAccessor.__getattr__(self, name)
    115 # Look up _geometry_type_name via __dict__ directly so that this method
    116 # cannot recurse when the instance has been constructed without
    117 # __init__ (e.g. during deepcopy / unpickling).
    118 geom_name = self.__dict__.get("_geometry_type_name", "<unknown>")
--> 119 raise AttributeError(
    120     f"DataArray with geometry '{geom_name}' has no z-coordinates; "
    121     "only layered 3D dfsu DataArrays expose .z"
    122 )

AttributeError: DataArray with geometry 'Grid2D' has no z-coordinates; only layered 3D dfsu DataArrays expose .z

Creating a dfsu from a mesh

import pandas as pd

msh = mikeio.open("../data/odense_rough.mesh")
data = 36.5 * np.ones(msh.geometry.n_elements)
time = pd.date_range("2020-1-1", periods=1)
item = mikeio.ItemInfo(mikeio.EUMType.Chezy_No)

da = mikeio.DataArray(data, time=time, item=item, geometry=msh.geometry)
da.to_dfs("chezy_map.dfsu")
import os
os.remove("chezy_map.dfsu")