Merge pull request #691 from DHI/append

Dfs append data to existing file

mikeio/dfs/_dfs2.py

@@ -242,6 +242,37 @@ def read(
             validate=False,
         )
 
+    def append(self, ds: Dataset, validate: bool = True) -> None:
+        """
+        Append a Dataset to an existing dfs2 file
+
+        Parameters
+        ----------
+        ds: Dataset
+            Dataset to append
+        validate: bool, optional
+            Check if the dataset to append has the same geometry and items as the original file,
+            by default True
+
+        Notes
+        -----
+        The original file is modified.
+        """
+        if validate:
+            if self.geometry != ds.geometry:
+                raise ValueError("The geometry of the dataset to append does not match")
+
+            for item_s, item_o in zip(ds.items, self.items):
+                if item_s != item_o:
+                    raise ValueError(
+                        f"Item in dataset {item_s.name} does not match {item_o.name}"
+                    )
+
+        dfs = DfsFileFactory.Dfs2FileOpenAppend(str(self._filename))
+        write_dfs_data(dfs=dfs, ds=ds, n_spatial_dims=2)
+
+        self._n_timesteps = dfs.FileInfo.TimeAxis.NumberOfTimeSteps
+
     def _open(self) -> None:
         self._dfs = DfsFileFactory.Dfs2FileOpen(self._filename)
         self._source = self._dfs

mikeio/dfs/_dfs3.py

@@ -267,6 +267,35 @@ def read(
             validate=False,
         )
 
+    def append(self, ds: Dataset, validate: bool = True) -> None:
+        """
+        Append a Dataset to an existing dfs3 file
+
+        Parameters
+        ----------
+        ds: Dataset
+            Dataset to append
+        validate: bool, optional
+            Validate that the dataset to append has the same geometry and items as the original file
+
+        Notes
+        -----
+        The original file is modified.
+        """
+        if validate:
+            if self.geometry != ds.geometry:
+                raise ValueError("The geometry of the dataset to append does not match")
+
+            for item_s, item_o in zip(ds.items, self.items):
+                if item_s != item_o:
+                    raise ValueError(
+                        f"Item in dataset {item_s.name} does not match {item_o.name}"
+                    )
+
+        dfs = DfsFileFactory.Dfs3FileOpenAppend(str(self._filename))
+        write_dfs_data(dfs=dfs, ds=ds, n_spatial_dims=3)
+        self._n_timesteps = dfs.FileInfo.TimeAxis.NumberOfTimeSteps
+
     @staticmethod
     def _get_bottom_values(data: np.ndarray) -> np.ndarray:
 

mikeio/dfsu/_dfsu.py

@@ -9,6 +9,7 @@
 from mikecore.DfsFactory import DfsFactory
 from mikecore.DfsuBuilder import DfsuBuilder
 from mikecore.DfsuFile import DfsuFile, DfsuFileType
+from mikecore.DfsFileFactory import DfsFileFactory
 from mikecore.eum import eumQuantity, eumUnit
 from tqdm import trange
 
@@ -47,8 +48,6 @@ def write_dfsu(filename: str | Path, data: Dataset) -> None:
         raise ValueError("Non-equidistant time axis is not supported.")
 
     dt = data.timestep
-    n_time_steps = len(data.time)
-
     geometry = data.geometry
     dfsu_filetype = DfsuFileType.Dfsu2D
 
@@ -84,8 +83,16 @@ def write_dfsu(filename: str | Path, data: Dataset) -> None:
     builder.ApplicationVersion = __dfs_version__
     dfs = builder.CreateFile(filename)
 
+    write_dfsu_data(dfs, data, geometry.is_layered)
+
+
+def write_dfsu_data(dfs: DfsuFile, ds: Dataset, is_layered: bool) -> None:
+
+    n_time_steps = len(ds.time)
+    data = ds
+
     for i in range(n_time_steps):
-        if geometry.is_layered:
+        if is_layered:
             if "time" in data.dims:
                 assert data._zn is not None
                 zn = data._zn[i]
@@ -492,12 +499,40 @@ def read(
             dt=self.timestep,
         )
 
+
+    def append(self, ds: Dataset, validate: bool = True) -> None:
+        """
+        Append data to an existing dfsu file
+
+        Parameters
+        ----------
+        data: Dataset
+            Dataset to be appended
+        validate: bool, optional
+            Validate that the items and geometry match, by default True
+        """
+        if validate:
+            if ds.geometry != self.geometry:
+                raise ValueError("The geometry of the dataset to append does not match")
+
+            for item_s, item_o in zip(ds.items, self.items):
+                if item_s != item_o:
+                    raise ValueError(
+                        f"Item in dataset {item_s.name} does not match {item_o.name}"
+                    )
+
+        dfs = DfsFileFactory.DfsuFileOpenAppend(str(self._filename), parameters=None)
+        write_dfsu_data(dfs=dfs, ds=ds, is_layered=False)
+        info = _get_dfsu_info(self._filename)
+        self._time = info.time
+
     def _parse_geometry_sel(
         self,
         area: Tuple[float, float, float, float] | None,
         x: float | None,
         y: float | None,
     ) -> np.ndarray | None:
+
         """Parse geometry selection
 
         Parameters

mikeio/dfsu/_layered.py

@@ -5,6 +5,7 @@
 from matplotlib.axes import Axes
 import numpy as np
 from mikecore.DfsuFile import DfsuFile, DfsuFileType
+from mikecore.DfsFileFactory import DfsFileFactory
 import pandas as pd
 from scipy.spatial import cKDTree
 from tqdm import trange
@@ -30,6 +31,7 @@
     get_nodes_from_source,
     get_elements_from_source,
     _validate_elements_and_geometry_sel,
+    write_dfsu_data,
 )
 
 if TYPE_CHECKING:
@@ -363,6 +365,32 @@ def read(
                 dt=self.timestep,
             )
 
+    def append(self, ds: Dataset, validate: bool = True) -> None:
+        """
+        Append data to a dfsu file
+
+        Parameters
+        ---------
+        ds: Dataset
+            Dataset to append
+        validate: bool, optional
+            Validate that the dataset to append has the same geometry and items, by default True
+        """
+        if validate:
+            if self.geometry != ds.geometry:
+                raise ValueError("The geometry of the dataset to append does not match")
+
+            for item_s, item_o in zip(ds.items, self.items):
+                if item_s != item_o:
+                    raise ValueError(
+                        f"Item in dataset {item_s.name} does not match {item_o.name}"
+                    )
+
+        dfs = DfsFileFactory.DfsuFileOpenAppend(str(self._filename), parameters=None)
+        write_dfsu_data(dfs=dfs, ds=ds, is_layered=ds.geometry.is_layered)
+        info = _get_dfsu_info(self._filename)
+        self._time = info.time
+
 
 class Dfsu2DV(DfsuLayered):
     def plot_vertical_profile(

mikeio/spatial/_FM_geometry.py

@@ -386,6 +386,21 @@ def boundary_codes(self) -> list[int]:
         valid = list(set(self.codes))
         return [code for code in valid if code > 0]
 
+    def __eq__(self, value: Any) -> bool:
+
+        # this is not an exhaustive check, but should be sufficient for most cases
+
+        if self.__class__ != value.__class__:
+            return False
+
+        if self.projection != value.projection:
+            return False
+
+        if not np.array_equal(self.node_coordinates, value.node_coordinates):
+            return False
+
+        return True
+
 
 class GeometryFM2D(_GeometryFM):
     def __init__(

tests/test_dataarray.py

@@ -1270,7 +1270,7 @@ def test_xzy_selection():
     das_xzy = ds.Temperature.sel(x=348946, y=6173673, z=0)
 
     # check for point geometry after selection
-    assert type(das_xzy.geometry) == mikeio.spatial.GeometryPoint3D
+    assert type(das_xzy.geometry) is mikeio.spatial.GeometryPoint3D
     assert das_xzy.values[0] == pytest.approx(17.381)
 
     # do the same but go one level deeper, but finding the index first
@@ -1301,7 +1301,7 @@ def test_layer_selection():
 
     das_layer = ds.Temperature.sel(layers=0)
     # should not be layered after selection
-    assert type(das_layer.geometry) == mikeio.spatial.GeometryFM2D
+    assert type(das_layer.geometry) is mikeio.spatial.GeometryFM2D
 
 
 def test_time_selection():

tests/test_dataset.py

@@ -1504,7 +1504,7 @@ def test_layer_selection():
 
     dss_layer = ds.sel(layers=0)
     # should not be layered after selection
-    assert type(dss_layer.geometry) == mikeio.spatial.GeometryFM2D
+    assert type(dss_layer.geometry) is mikeio.spatial.GeometryFM2D
 
 
 def test_time_selection():

tests/test_dfs2.py

@@ -864,3 +864,45 @@ def test_to_xarray():
     xr_da = da.to_xarray()
     assert xr_da.x[0] == pytest.approx(0.25)
     assert xr_da.y[0] == pytest.approx(0.25)
+
+
+def test_append(tmp_path):
+    ds = mikeio.read("tests/testdata/eq.dfs2", time=[0, 1])
+    ds2 = mikeio.read("tests/testdata/eq.dfs2", time=[2, 3])
+
+    new_filename = tmp_path / "eq_appended.dfs2"
+    ds.to_dfs(new_filename)
+    dfs = mikeio.open(new_filename)
+    assert dfs.n_timesteps == 2
+    dfs.append(ds2)
+    assert dfs.n_timesteps == 4
+
+    ds3 = mikeio.read(new_filename)
+    assert ds3.n_timesteps == 4
+    assert ds3.time[-1] == ds2.time[-1]
+    assert ds3[0].values[-1, 0, 0] == ds2[0].values[-1, 0, 0]
+
+
+def test_append_mismatch_items_not_possible(tmp_path):
+    ds = mikeio.read("tests/testdata/consistency/oresundHD.dfs2", time=[0, 1])
+    ds2 = mikeio.read(
+        "tests/testdata/consistency/oresundHD.dfs2", time=[2, 3], items=[1, 2]
+    )
+
+    new_filename = tmp_path / "eq_appended.dfs2"
+    ds.to_dfs(new_filename)
+    dfs = mikeio.open(new_filename)
+    with pytest.raises(ValueError, match="Item"):
+        dfs.append(ds2)
+
+
+def test_append_mismatch_geometry(tmp_path):
+    ds = mikeio.read("tests/testdata/consistency/oresundHD.dfs2", time=[0, 1])
+    ds2 = mikeio.read("tests/testdata/consistency/oresundHD.dfs2").isel(
+        x=slice(1, None)
+    )
+    new_filename = tmp_path / "append_mismatch_geometry.dfs2"
+    ds.to_dfs(new_filename)
+    dfs = mikeio.open(new_filename)
+    with pytest.raises(ValueError, match="geometry"):
+        dfs.append(ds2)

tests/test_dfs3.py

@@ -266,3 +266,16 @@ def test_to_xarray() -> None:
     xr_da = da.to_xarray()
     assert xr_da.x[0] == pytest.approx(0.25)
     assert xr_da.y[0] == pytest.approx(0.25)
+
+
+def test_append_dfs3(tmp_path):
+    fn = "tests/testdata/Karup_MIKE_SHE_head_output.dfs3"
+    ds = mikeio.read(fn, time=[0, 1])
+    new_fp = tmp_path / "test_append.dfs3"
+    ds.to_dfs(new_fp)
+
+    ds2 = mikeio.read(fn, time=[2, 3])
+
+    dfs = mikeio.open(new_fp)
+
+    dfs.append(ds2)

tests/test_dfsu.py

@@ -970,3 +970,21 @@ def test_writing_non_equdistant_dfsu_is_not_possible(tmp_path):
     with pytest.raises(ValueError, match="equidistant"):
         fp = tmp_path / "not_gonna_work.dfsu"
         dss.to_dfs(fp)
+
+
+def test_append_dfsu_2d(tmp_path):
+    ds = mikeio.read("tests/testdata/consistency/oresundHD.dfsu", time=[0, 1])
+    ds2 = mikeio.read("tests/testdata/consistency/oresundHD.dfsu", time=[2, 3])
+    new_filename = tmp_path / "appended.dfsu"
+    ds.to_dfs(new_filename)
+    dfs = mikeio.open(new_filename)
+    assert dfs.time[-1] == ds.time[-1]
+    dfs.append(ds2)
+    assert dfs.time[-1] == ds2.time[-1]
+
+    ds3 = mikeio.read(new_filename)
+    assert ds3.n_timesteps == 4
+    assert ds3.time[-1] == ds2.time[-1]
+    assert (
+        ds3.V_velocity.isel(time=3).values[0] == ds2.V_velocity.isel(time=1).values[0]
+    )

tests/test_dfsu_layered.py

@@ -617,6 +617,24 @@ def test_read_wildcard_items():
     assert ds.n_items == 1
 
 
+def test_append_dfsu_3d(tmp_path):
+    ds = mikeio.read("tests/testdata/oresund_sigma_z.dfsu", time=[0])
+    assert ds.timestep == pytest.approx(10800)
+    ds2 = mikeio.read("tests/testdata/oresund_sigma_z.dfsu", time=[1])
+    new_filename = tmp_path / "appended.dfsu"
+    ds.to_dfs(new_filename)
+    dfs = mikeio.open(new_filename)
+    assert dfs.timestep == pytest.approx(10800)
+    assert dfs.n_timesteps == 1
+    dfs.append(ds2)
+    assert dfs.n_timesteps == 2
+    assert dfs.time[-1] == ds2.time[-1]
+
+    # verify that the new file can be read
+    ds3 = mikeio.read(new_filename)
+    assert ds3.n_timesteps == 2
+    assert ds3.time[-1] == ds2.time[-1]
+
 def test_read_elements_3d():
     ds = mikeio.read("tests/testdata/oresund_sigma_z.dfsu", elements=[0, 10])
     assert ds.geometry.element_coordinates[0][0] == pytest.approx(354020.46382194717)
@@ -625,3 +643,4 @@ def test_read_elements_3d():
     ds2 = mikeio.read("tests/testdata/oresund_sigma_z.dfsu", elements=[10, 0])
     assert ds2.geometry.element_coordinates[1][0] == pytest.approx(354020.46382194717)
     assert ds2.Salinity.to_numpy()[0, 1] == pytest.approx(23.18906021118164)
+