Merge pull request #728 from DHI/directions-conversion

Fix bug in Dfsu spectral files directions conversion

mikeio/dfsu/_spectral.py

@@ -5,10 +5,11 @@
 import numpy as np
 import pandas as pd
 from mikecore.DfsuFile import DfsuFile, DfsuFileType
+from mikecore.DfsFileFactory import DfsFileFactory
 from tqdm import trange
 
 from ..dataset import DataArray, Dataset
-from ..eum import ItemInfo
+from ..eum import ItemInfo, EUMUnit
 from ..dfs._dfs import _get_item_info, _valid_item_numbers, _valid_timesteps
 from .._spectral import calc_m0_from_spectrum
 from ._dfsu import (
@@ -134,8 +135,12 @@ def _read_geometry(
         dfs = DfsuFile.Open(filename)
         dfsu_type = DfsuFileType(dfs.DfsuFileType)
 
-        dir = dfs.Directions
-        directions = None if dir is None else dir * (180 / np.pi)
+        directions = dfs.Directions
+        if directions is not None:
+            dir_unit = DfsuSpectral._get_direction_unit(filename)
+            dir_conversion = 180.0 / np.pi if dir_unit == int(EUMUnit.radian) else 1.0
+            directions = directions * dir_conversion
+
         frequencies = dfs.Frequencies
 
         # geometry
@@ -177,6 +182,19 @@ def _read_geometry(
         dfs.Close()
         return geometry
 
+    @staticmethod
+    def _get_direction_unit(filename: str) -> int:
+        """Determine if the directional axis is in degrees or radians"""
+        source = DfsFileFactory.DfsGenericOpen(filename)
+        try:
+            for static_item in iter(source.ReadStaticItemNext, None):
+                if static_item.Name == "Direction":
+                    return static_item.Quantity.Unit.value
+        finally:
+            source.Close()
+
+        raise ValueError("Direction static item not found in the file.")
+
     @property
     def n_frequencies(self) -> int | None:
         """Number of frequencies"""
@@ -279,15 +297,15 @@ def read(
 
         Examples
         --------
-        >>> mikeio.read("tests/testdata/line_spectra.dfsu")
+        >>> mikeio.read("tests/testdata/spectra/line_spectra.dfsu")
         <mikeio.Dataset>
         dims: (time:4, node:10, direction:16, frequency:25)
         time: 2017-10-27 00:00:00 - 2017-10-27 05:00:00 (4 records)
         geometry: DfsuSpectral1D (9 elements, 10 nodes)
         items:
           0:  Energy density <Wave energy density> (meter pow 2 sec per deg)
 
-        >>> mikeio.read("tests/testdata/area_spectra.dfsu", time=-1)
+        >>> mikeio.read("tests/testdata/spectra/area_spectra.dfsu", time=-1)
         <mikeio.Dataset>
         dims: (element:40, direction:16, frequency:25)
         time: 2017-10-27 05:00:00 (time-invariant)

notebooks/Dfs2 - Various types.ipynb

@@ -599,7 +599,7 @@
     }
    ],
    "source": [
-    "fn = \"../tests/testdata/dir_wave_analysis_spectra.dfs2\"\n",
+    "fn = \"../tests/testdata/spectra/dir_wave_analysis_spectra.dfs2\"\n",
     "dfs = mikeio.open(fn, type=\"spectral\")\n",
     "da = dfs.read()[0]\n",
     "da"
@@ -676,7 +676,7 @@
     }
    ],
    "source": [
-    "fn = \"../tests/testdata/pt_spectra.dfs2\"\n",
+    "fn = \"../tests/testdata/spectra/pt_spectra.dfs2\"\n",
     "dfs = mikeio.open(fn, type=\"spectral\")\n",
     "da = dfs.read()[0]\n",
     "da"

tests/test_dfs2.py

@@ -27,13 +27,13 @@ def dfs2_random_2items():
 
 @pytest.fixture
 def dfs2_pt_spectrum():
-    filepath = Path("tests/testdata/pt_spectra.dfs2")
+    filepath = Path("tests/testdata/spectra/pt_spectra.dfs2")
     return mikeio.open(filepath, type="spectral")
 
 
 @pytest.fixture
 def dfs2_pt_spectrum_linearf():
-    filepath = Path("tests/testdata/dir_wave_analysis_spectra.dfs2")
+    filepath = Path("tests/testdata/spectra/dir_wave_analysis_spectra.dfs2")
     return mikeio.open(filepath, type="spectral")
 
 
@@ -58,7 +58,6 @@ def test_get_time_without_reading_data():
 
 
 def test_write_projected(tmp_path):
-
     fp = tmp_path / "utm.dfs2"
 
     nt = 100
@@ -161,7 +160,6 @@ def test_write_without_time(tmp_path):
 
 
 def test_read(dfs2_random):
-
     dfs = dfs2_random
     assert isinstance(dfs.geometry, Grid2D)
     ds = dfs.read(items=["testing water level"])
@@ -180,7 +178,6 @@ def test_read_bad_item(dfs2_random):
 
 
 def test_read_temporal_subset_slice():
-
     filename = r"tests/testdata/eq.dfs2"
     dfs = mikeio.open(filename)
     ds = dfs.read(time=slice("2000-01-01 00:00", "2000-01-01 12:00"))
@@ -189,7 +186,6 @@ def test_read_temporal_subset_slice():
 
 
 def test_read_area_subset_bad_bbox():
-
     filename = "tests/testdata/europe_wind_long_lat.dfs2"
     bbox = (10, 40, 20)
     with pytest.raises(ValueError):
@@ -220,7 +216,6 @@ def test_subset_bbox():
 
 
 def test_read_area_subset():
-
     filename = "tests/testdata/eq.dfs2"
     bbox = [10, 4, 12, 7]
 
@@ -247,7 +242,6 @@ def test_read_area_subset():
 
 
 def test_read_numbered_access(dfs2_random_2items):
-
     dfs = dfs2_random_2items
 
     res = dfs.read(items=[1])
@@ -349,7 +343,7 @@ def test_properties_pt_spectrum_linearf(dfs2_pt_spectrum_linearf):
 
 def test_dir_wave_spectra_relative_time_axis():
     ds = mikeio.open(
-        "tests/testdata/dir_wave_analysis_spectra.dfs2", type="spectral"
+        "tests/testdata/spectra/dir_wave_analysis_spectra.dfs2", type="spectral"
     ).read()
     assert ds.n_items == 1
     assert ds.geometry.nx == 128
@@ -431,7 +425,6 @@ def test_select_area_rotated_UTM_2():
 
 
 def test_write_selected_item_to_new_file(dfs2_random_2items, tmp_path):
-
     dfs = dfs2_random_2items
 
     fp = tmp_path / "simple.dfs2"
@@ -455,7 +448,6 @@ def test_write_selected_item_to_new_file(dfs2_random_2items, tmp_path):
 
 
 def test_repr(dfs2_gebco):
-
     text = repr(dfs2_gebco)
 
     assert "Dfs2" in text
@@ -464,7 +456,6 @@ def test_repr(dfs2_gebco):
 
 
 def test_repr_time(dfs2_random):
-
     dfs = dfs2_random
     text = repr(dfs)
 
@@ -475,7 +466,6 @@ def test_repr_time(dfs2_random):
 
 
 def test_write_modified_data_to_new_file(dfs2_gebco, tmp_path):
-
     dfs = dfs2_gebco
 
     fp = tmp_path / "mod.dfs2"
@@ -492,7 +482,6 @@ def test_write_modified_data_to_new_file(dfs2_gebco, tmp_path):
 
 
 def test_read_some_time_step(dfs2_random_2items):
-
     dfs = dfs2_random_2items
     res = dfs.read(time=[1, 2])
 
@@ -501,7 +490,6 @@ def test_read_some_time_step(dfs2_random_2items):
 
 
 def test_interpolate_non_equidistant_data(tmp_path):
-
     ds = mikeio.read(
         "tests/testdata/eq.dfs2", time=[0, 2, 3, 6]
     )  # non-equidistant dataset
@@ -525,7 +513,6 @@ def test_interpolate_non_equidistant_data(tmp_path):
 
 
 def test_write_some_time_step(tmp_path):
-
     ds = mikeio.read("tests/testdata/waves.dfs2", time=[1, 2])
 
     assert ds[0].to_numpy().shape[0] == 2
@@ -593,7 +580,6 @@ def test_write_accumulated_datatype(tmp_path):
 
 
 def test_write_NonEqCalendarAxis(tmp_path):
-
     fp = tmp_path / "simple.dfs2"
 
     d = np.random.random([6, 5, 10])
@@ -623,7 +609,6 @@ def test_write_NonEqCalendarAxis(tmp_path):
 
 
 def test_write_non_equidistant_data(tmp_path):
-
     ds = mikeio.read(
         "tests/testdata/eq.dfs2", time=[0, 2, 3, 6]
     )  # non-equidistant dataset
@@ -658,7 +643,6 @@ def test_read_concat_write_dfs2(tmp_path):
 
 
 def test_spatial_aggregation_dfs2_to_dfs0(tmp_path):
-
     outfilename = tmp_path / "waves_max.dfs0"
 
     ds = mikeio.read("tests/testdata/waves.dfs2")
@@ -705,15 +689,13 @@ def test_grid2d_plot():
 
 
 def test_read_single_precision():
-
     ds = mikeio.read("tests/testdata/random.dfs2", items=0, dtype=np.float32)
 
     assert len(ds) == 1
     assert ds[0].dtype == np.float32
 
 
 def dfs2_props_to_list(d):
-
     lon = d._dfs.FileInfo.Projection.Longitude
     lat = d._dfs.FileInfo.Projection.Latitude
     rot = d._dfs.FileInfo.Projection.Orientation
@@ -788,7 +770,7 @@ def test_read_write_header_unchanged_vertical(tmp_path):
 
 
 def test_read_write_header_unchanged_spectral_2(tmp_path):
-    is_header_unchanged_on_read_write(tmp_path, "pt_spectra.dfs2")
+    is_header_unchanged_on_read_write(tmp_path, "spectra/pt_spectra.dfs2")
 
 
 def test_read_write_header_unchanged_MIKE_SHE_output(tmp_path):
@@ -812,7 +794,6 @@ def test_MIKE_SHE_output():
 
 
 def test_read_dfs2_static_dt_zero():
-
     with pytest.warns(UserWarning, match="positive"):
         ds = mikeio.read("tests/testdata/single_time_dt_zero.dfs2")
     assert ds.n_timesteps == 1

tests/test_dfsu_plot.py

@@ -135,7 +135,6 @@ def test_da_plot():
 
 
 def test_plot_non_utm_file():
-
     ds = mikeio.read("tests/testdata/FakeLake_NONUTM.dfsu")
     da = ds[0]
     da.plot()
@@ -159,16 +158,16 @@ def test_plot_vertical_transect():
 
 def test_plot_point_spectrum():
     # directional spectra
-    da = mikeio.read("tests/testdata/line_dir_spectra.dfsu")[0]
+    da = mikeio.read("tests/testdata/spectra/line_dir_spectra.dfsu")[0]
     da.isel(node=4).plot()
 
     # frequency spectra
-    da2 = mikeio.read("tests/testdata/line_freq_spectra.dfsu")[0]
+    da2 = mikeio.read("tests/testdata/spectra/line_freq_spectra.dfsu")[0]
     da2_pt = da2.isel(node=4)
     da2_pt.plot()
 
     # 2d spectra
-    da_pt = mikeio.read("tests/testdata/pt_spectra.dfsu")[0]
+    da_pt = mikeio.read("tests/testdata/spectra/pt_spectra.dfsu")[0]
     da_pt.plot.patch()
     da_pt.plot.contour()
     da_pt.plot.contourf()

tests/test_dfsu_spectral.py

@@ -10,43 +10,49 @@
 
 @pytest.fixture
 def dfsu_pt():
-    filename = "tests/testdata/pt_spectra.dfsu"
+    filename = "tests/testdata/spectra/pt_spectra.dfsu"
     return mikeio.open(filename)
 
 
 @pytest.fixture
 def dfsu_line():
-    filename = "tests/testdata/line_spectra.dfsu"
+    filename = "tests/testdata/spectra/line_spectra.dfsu"
+    return mikeio.open(filename)
+
+
+@pytest.fixture
+def dfsu_line_degrees():
+    filename = "tests/testdata/spectra/line_spectra_degrees.dfsu"
     return mikeio.open(filename)
 
 
 @pytest.fixture
 def dfsu_area():
-    filename = "tests/testdata/area_spectra.dfsu"
+    filename = "tests/testdata/spectra/area_spectra.dfsu"
     return mikeio.open(filename)
 
 
 @pytest.fixture
 def dfsu_area_sector():
-    filename = "tests/testdata/MIKE21SW_dir_sector_area_spectra.dfsu"
+    filename = "tests/testdata/spectra/MIKE21SW_dir_sector_area_spectra.dfsu"
     return mikeio.open(filename)
 
 
 @pytest.fixture
 def dfsu_pt_freq():
-    filename = "tests/testdata/pt_freq_spectra.dfsu"
+    filename = "tests/testdata/spectra/pt_freq_spectra.dfsu"
     return mikeio.open(filename)
 
 
 @pytest.fixture
 def dfsu_line_dir():
-    filename = "tests/testdata/line_dir_spectra.dfsu"
+    filename = "tests/testdata/spectra/line_dir_spectra.dfsu"
     return mikeio.open(filename)
 
 
 @pytest.fixture
 def dfsu_area_freq():
-    filename = "tests/testdata/area_freq_spectra.dfsu"
+    filename = "tests/testdata/spectra/area_freq_spectra.dfsu"
     return mikeio.open(filename)
 
 
@@ -70,6 +76,18 @@ def test_properties_line_spectrum(dfsu_line):
     assert len(dfs.directions) == 16
     assert dfs.geometry.n_nodes == 10
     assert dfs.geometry.n_elements == 9
+    dir = dfs.geometry.directions
+    assert dir[0] == pytest.approx(0.0)
+    assert dir[-1] == pytest.approx(337.5)
+
+
+def test_properties_line_spectrum_degrees(dfsu_line_degrees):
+    dfs = dfsu_line_degrees
+    assert dfs.geometry.is_spectral
+    assert dfs._type == DfsuFileType.DfsuSpectral1D
+    dir = dfs.geometry.directions
+    assert dir[0] == pytest.approx(0.0)
+    assert dir[-1] == pytest.approx(350.0)
 
 
 def test_properties_area_spectrum(dfsu_area):