Merge pull request #837 from DHI/feature/write_spectral_dfsu

Feature/write_spectral_dfsu

Author: ecomodeller

mikeio/dfsu/_dfsu.py

@@ -51,7 +51,7 @@ def write_dfsu(filename: str | Path, data: Dataset) -> None:
     geometry = data.geometry
     dfsu_filetype = DfsuFileType.Dfsu2D
 
-    if geometry.is_layered:
+    if geometry.is_layered or geometry.is_spectral:
         dfsu_filetype = geometry._type.value
 
     xn = geometry.node_coordinates[:, 0]
@@ -61,9 +61,24 @@ def write_dfsu(filename: str | Path, data: Dataset) -> None:
     elem_table = [np.array(e) + 1 for e in geometry.element_table]
 
     builder = DfsuBuilder.Create(dfsu_filetype)
-    if dfsu_filetype != DfsuFileType.Dfsu2D:
+    if dfsu_filetype in (
+        DfsuFileType.Dfsu3DSigma,
+        DfsuFileType.Dfsu3DSigmaZ,
+        DfsuFileType.DfsuVerticalColumn,
+        DfsuFileType.DfsuVerticalProfileSigma,
+        DfsuFileType.DfsuVerticalProfileSigmaZ,
+    ):
         builder.SetNumberOfSigmaLayers(geometry.n_sigma_layers)
 
+    if dfsu_filetype in (
+        DfsuFileType.DfsuSpectral0D,
+        DfsuFileType.DfsuSpectral1D,
+        DfsuFileType.DfsuSpectral2D,
+    ):
+        builder.SetFrequencies(geometry.frequencies)
+        # TODO should directions always be converted to radians
+        builder.SetDirections(np.deg2rad(geometry.directions))
+
     builder.SetNodes(xn, yn, zn, geometry.codes)
     builder.SetElements(elem_table)
 
@@ -82,9 +97,6 @@ def write_dfsu(filename: str | Path, data: Dataset) -> None:
     builder.SetTemporalAxis(temporal_axis)
     builder.SetZUnit(eumUnit.eumUmeter)
 
-    if dfsu_filetype != DfsuFileType.Dfsu2D:
-        builder.SetNumberOfSigmaLayers(geometry.n_sigma_layers)
-
     for item in data.items:
         builder.AddDynamicItem(item.name, eumQuantity.Create(item.type, item.unit))
 
@@ -117,7 +129,13 @@ def write_dfsu_data(dfs: DfsuFile, ds: Dataset, is_layered: bool) -> None:
                 d = da.to_numpy()[i, :]
             else:
                 d = da.to_numpy()
+            d = d.copy()  # to avoid modifying the input
             d[np.isnan(d)] = data.deletevalue
+
+            if d.ndim != 1:
+                d = np.moveaxis(d, 0, -1)
+                d = d.reshape(-1)
+
             dfs.WriteItemTimeStepNext(t_rel[i], d.astype(np.float32))
     dfs.Close()
 

mikeio/spatial/_FM_geometry_spectral.py

@@ -105,6 +105,10 @@ def n_frequencies(self) -> int:
         """Number of frequencies."""
         return 0 if self.frequencies is None else len(self.frequencies)
 
+    @property
+    def is_layered(self) -> bool:
+        return False
+
     @property
     def frequencies(self) -> np.ndarray | None:
         """Frequency axis."""
@@ -182,6 +186,31 @@ def elements_to_geometry(  # type: ignore
 class GeometryFMLineSpectrum(_GeometryFMSpectrum):
     """Flexible mesh line spectrum geometry."""
 
+    @staticmethod
+    def create_dummy_coordinates(
+        n_nodes: int, frequencies: np.ndarray, directions: np.ndarray
+    ) -> GeometryFMLineSpectrum:
+        # boogus x, y, z 1..n, 1..n, -10
+        node_coordinates = np.zeros((n_nodes, 3), dtype=float)
+        node_coordinates[:, 0] = np.arange(1, n_nodes + 1)  # x
+        node_coordinates[:, 1] = np.arange(1, n_nodes + 1)  # y
+        node_coordinates[:, 2] = -10.0  # z
+
+        # element table with nodes [0,1], [1,2], ..., [n-2,n-1]
+        element_table = np.zeros((n_nodes - 1, 2), dtype=int)
+        for i in range(n_nodes - 1):
+            element_table[i, 0] = i
+            element_table[i, 1] = i + 1
+        return GeometryFMLineSpectrum(
+            node_coordinates=node_coordinates,
+            element_table=element_table,
+            codes=np.zeros(n_nodes, dtype=int),
+            dfsu_type=DfsuFileType.DfsuSpectral1D,
+            projection="LONG/LAT",
+            frequencies=frequencies,
+            directions=directions,
+        )
+
     def isel(  # type: ignore
         self, idx: Sequence[int], axis: str = "node"
     ) -> GeometryFMPointSpectrum | GeometryFMLineSpectrum:

pyproject.toml

@@ -20,9 +20,9 @@ exclude = [
 
 [project]
 name = "mikeio"
-version = "2.7.0"
+version = "3.0.0a0"
 dependencies = [
-    "mikecore>=0.2.1",
+    "mikecore>=0.3.0a0",
     "numpy>=1.22.0",
     "pandas>=1.3",
     "matplotlib>=3.6.0",

tests/test_dfsu_spectral.py

@@ -1,12 +1,14 @@
+from pathlib import Path
 import pytest
 import numpy as np
 import matplotlib.pyplot as plt
 import mikeio
 from mikecore.DfsuFile import DfsuFileType
 
-from mikeio import DfsuSpectral
+from mikeio import DfsuSpectral, EUMType, EUMUnit
 from mikeio.spatial import GeometryFMPointSpectrum, GeometryFMAreaSpectrum
 import mikeio._spectral as _spectral
+from mikeio.spatial._FM_geometry_spectral import GeometryFMLineSpectrum
 
 
 @pytest.fixture
@@ -439,3 +441,67 @@ def test_plot_da_spectrum(dfsu_pt: DfsuSpectral) -> None:
     # dfs.plot_spectrum(spec, title="pt", plot_type="shaded")
     # dfs.plot_spectrum(spec, r_as_periods=False, plot_type="contour")
     plt.close("all")
+
+
+def test_write_line_spectra(dfsu_line: DfsuSpectral, tmp_path: Path) -> None:
+    ds = dfsu_line.read()
+
+    fp = tmp_path / "line.dfsu"
+
+    ds.to_dfs(fp)
+
+    ds2 = mikeio.read(fp)
+
+
+def test_write_line_spectra_energy(tmp_path: Path) -> None:
+    ds = mikeio.read("tests/testdata/spectra/North_BC_2024_subset.dfsu")
+    assert ds["Energy density"].type == EUMType.Wave_energy_density
+    assert ds["Energy density"].unit == EUMUnit.meter_pow_2_sec_per_deg
+    assert ds.geometry.directions[-1] == pytest.approx(350)
+    assert ds.geometry.frequencies[0] == pytest.approx(0.035)
+    assert ds.geometry.frequencies[-1] == pytest.approx(0.983586)
+    assert ds["Energy density"].isel(time=-1).isel(node=0).isel(frequency=35).isel(
+        direction=28
+    ).to_numpy() == pytest.approx(2.6239042e-05)
+
+    fp = tmp_path / "line_energy.dfsu"
+    ds.to_dfs(fp)
+
+    ds2 = mikeio.read(fp)
+    assert ds2["Energy density"].type == EUMType.Wave_energy_density
+    assert ds2["Energy density"].unit == EUMUnit.meter_pow_2_sec_per_deg
+    assert ds2.geometry.directions[-1] == pytest.approx(350)
+    assert ds2.geometry.frequencies[0] == pytest.approx(0.035)
+    assert ds2.geometry.frequencies[-1] == pytest.approx(0.983586)
+    assert ds2["Energy density"].isel(time=-1).isel(node=0).isel(frequency=35).isel(
+        direction=28
+    ).to_numpy() == pytest.approx(2.6239042e-05)
+    assert np.all(ds2.to_numpy() == ds.to_numpy())
+
+
+def test_write_area_spectra(dfsu_area: DfsuSpectral, tmp_path: Path) -> None:
+    ds = dfsu_area.read()
+    fp = tmp_path / "area.dfsu"
+    ds.to_dfs(fp)
+
+    dfs = mikeio.DfsuSpectral(fp)
+    assert dfs.geometry.is_spectral
+    assert dfs._type == DfsuFileType.DfsuSpectral2D
+    assert dfs.geometry.n_frequencies == 25
+    assert dfs.frequencies is not None
+    assert len(dfs.frequencies) == 25
+    assert dfs.geometry.n_directions == 16
+    assert dfs.directions is not None
+    assert len(dfs.directions) == 16
+
+    ds2 = dfs.read()
+    assert np.all(ds2.to_numpy() == ds.to_numpy())
+
+
+def test_create_line_spectrum_dummy_coordinates() -> None:
+    freq = np.array([0.035, 0.983586])
+    dirs = np.arange(5.0, 360.0, step=10.0)
+    geometry = GeometryFMLineSpectrum.create_dummy_coordinates(
+        n_nodes=10, frequencies=freq, directions=dirs
+    )
+    assert geometry.n_nodes == 10