Repository: DHI/mikeio
Pull request #691 (2024-05-08 06:02:01)

Dfs append data to existing file
Common problem file ending on xxx-12-31 23:00

Let's extend the file with the last two timestep repeated, without the need of reading and writing the entire file.

![image](https://github.com/DHI/mikeio/assets/614215/33dc59f7-b077-4243-bfdc-a2909287f4b0)
![image](https://github.com/DHI/mikeio/assets/614215/e9a86b50-9321-4657-b10e-ebf9d5652c30)



Patch:
diff --git a/mikeio/dfs/_dfs2.py b/mikeio/dfs/_dfs2.py
index 6c5d8743..f656c3ec 100644
--- a/mikeio/dfs/_dfs2.py
+++ b/mikeio/dfs/_dfs2.py
@@ -242,6 +242,37 @@ class Dfs2(_Dfs123):
             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
diff --git a/mikeio/dfs/_dfs3.py b/mikeio/dfs/_dfs3.py
index 04e8f84a..9db45dfc 100644
--- a/mikeio/dfs/_dfs3.py
+++ b/mikeio/dfs/_dfs3.py
@@ -267,6 +267,35 @@ class Dfs3(_Dfs123):
             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:
 
diff --git a/mikeio/dfsu/_dfsu.py b/mikeio/dfsu/_dfsu.py
index 5e4f1e77..e8830227 100644
--- a/mikeio/dfsu/_dfsu.py
+++ b/mikeio/dfsu/_dfsu.py
@@ -9,6 +9,7 @@ import pandas as pd
 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 @@ class Dfsu2DH:
             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
diff --git a/mikeio/dfsu/_layered.py b/mikeio/dfsu/_layered.py
index 2519913c..9e0e46a7 100644
--- a/mikeio/dfsu/_layered.py
+++ b/mikeio/dfsu/_layered.py
@@ -5,6 +5,7 @@ from typing import Any, Sequence, Tuple, TYPE_CHECKING
 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 @@ from ._dfsu import (
     get_nodes_from_source,
     get_elements_from_source,
     _validate_elements_and_geometry_sel,
+    write_dfsu_data,
 )
 
 if TYPE_CHECKING:
@@ -363,6 +365,32 @@ class DfsuLayered:
                 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(
diff --git a/mikeio/spatial/_FM_geometry.py b/mikeio/spatial/_FM_geometry.py
index c930f3cf..cc56e905 100644
--- a/mikeio/spatial/_FM_geometry.py
+++ b/mikeio/spatial/_FM_geometry.py
@@ -386,6 +386,21 @@ class _GeometryFM(_Geometry):
         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__(


Test patch:
diff --git a/tests/test_dataarray.py b/tests/test_dataarray.py
index b40dc4c9..b1a81cba 100644
--- a/tests/test_dataarray.py
+++ b/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():
diff --git a/tests/test_dataset.py b/tests/test_dataset.py
index 82576fce..a5ac411a 100644
--- a/tests/test_dataset.py
+++ b/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():
diff --git a/tests/test_dfs2.py b/tests/test_dfs2.py
index e49adfae..39989858 100644
--- a/tests/test_dfs2.py
+++ b/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)
diff --git a/tests/test_dfs3.py b/tests/test_dfs3.py
index 9d7c2146..f941ce4f 100644
--- a/tests/test_dfs3.py
+++ b/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)
diff --git a/tests/test_dfsu.py b/tests/test_dfsu.py
index f3092f90..4d04a81a 100644
--- a/tests/test_dfsu.py
+++ b/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]
+    )
diff --git a/tests/test_dfsu_layered.py b/tests/test_dfsu_layered.py
index 5eaf93c1..d21dc957 100644
--- a/tests/test_dfsu_layered.py
+++ b/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)
+
diff --git a/tests/test_dfsu_transect.py b/tests/test_dfsu_transect.py
index 97d7ec4a..d4670fad 100644
--- a/tests/test_dfsu_transect.py
+++ b/tests/test_dfsu_transect.py
@@ -63,7 +63,7 @@ def test_transect_geometry_properties_geo(vslice_geo):
 def test_transect_read(vslice):
     ds = vslice.read()
     assert ds.geometry._type == DfsuFileType.DfsuVerticalProfileSigmaZ
-    assert type(ds.geometry) == GeometryFMVerticalProfile
+    assert type(ds.geometry) is GeometryFMVerticalProfile
     assert ds.n_items == 2
     assert ds.Salinity.name == "Salinity"
     assert ds.geometry.n_elements == 441
@@ -78,7 +78,7 @@ def test_transect_getitem_element(vslice):
     da = vslice.read().Salinity
     idx = 5
     da2 = da[:, idx]
-    assert type(da2.geometry) == GeometryPoint3D
+    assert type(da2.geometry) is GeometryPoint3D
     assert da2.geometry.x == da.geometry.element_coordinates[idx, 0]
     assert da2.geometry.y == da.geometry.element_coordinates[idx, 1]
     assert da2.geometry.z == da.geometry.element_coordinates[idx, 2]
@@ -90,7 +90,7 @@ def test_transect_isel(vslice):
     da = vslice.read().Salinity
     idx = 5
     da2 = da.isel(element=idx)
-    assert type(da2.geometry) == GeometryPoint3D
+    assert type(da2.geometry) is GeometryPoint3D
     assert da2.geometry.x == da.geometry.element_coordinates[idx, 0]
     assert da2.geometry.y == da.geometry.element_coordinates[idx, 1]
     assert da2.geometry.z == da.geometry.element_coordinates[idx, 2]
@@ -105,7 +105,7 @@ def test_transect_isel_multiple(vslice_geo):
 
     ds2 = ds.isel(element=idx)
     assert ds2.geometry.n_elements == 579
-    assert type(ds2.geometry) == GeometryFMVerticalProfile
+    assert type(ds2.geometry) is GeometryFMVerticalProfile
     rd2 = ds2.geometry.relative_element_distance
     assert rd2.max() < 15000
 
@@ -114,7 +114,7 @@ def test_transect_sel_time(vslice):
     da = vslice.read()[0]
     idx = 1
     da2 = da.sel(time="1997-09-16 00:00")
-    assert type(da2.geometry) == GeometryFMVerticalProfile
+    assert type(da2.geometry) is GeometryFMVerticalProfile
     assert np.all(da2.to_numpy() == da.to_numpy()[idx, :])
     assert da2.time[0] == da.time[1]
     assert da2.dims == ("element",)
@@ -124,7 +124,7 @@ def test_transect_sel_time(vslice):
 def test_transect_sel_xyz(vslice_geo):
     da = vslice_geo.read().Temperature
     da2 = da.sel(x=10.8, y=55.6, z=-3)
-    assert type(da2.geometry) == GeometryPoint3D
+    assert type(da2.geometry) is GeometryPoint3D
     assert da2.geometry.x == pytest.approx(10.802878)
     assert da2.geometry.y == pytest.approx(55.603096)
     assert da2.geometry.z == pytest.approx(-2.3)
@@ -134,13 +134,13 @@ def test_transect_sel_xyz(vslice_geo):
 def test_transect_sel_layers(vslice_geo):
     ds = vslice_geo.read()
     ds2 = ds.sel(layers=range(-6, -1))
-    assert type(ds2.geometry) == GeometryFMVerticalProfile
+    assert type(ds2.geometry) is GeometryFMVerticalProfile
 
 
 def test_transect_sel_xy(vslice_geo):
     da = vslice_geo.read().Temperature
     da2 = da.sel(x=10.8, y=55.6)
-    assert type(da2.geometry) == GeometryFMVerticalColumn
+    assert type(da2.geometry) is GeometryFMVerticalColumn
     gx, gy, _ = da2.geometry.element_coordinates.mean(axis=0)
     assert gx == pytest.approx(10.802878)
     assert gy == pytest.approx(55.603096)
diff --git a/tests/test_geometry_fm.py b/tests/test_geometry_fm.py
index 02ed5554..51425714 100644
--- a/tests/test_geometry_fm.py
+++ b/tests/test_geometry_fm.py
@@ -236,3 +236,44 @@ def test_layered(simple_3d_geom: GeometryFM3D):
 
     assert "elements: 2" in repr(g2)
     assert "layers: 2" in repr(g2)
+
+
+def test_equality():
+    nc = [
+        (0.0, 0.0, 0.0),  # 0
+        (1.0, 0.0, 0.0),  # 1
+        (1.0, 1.0, 0.0),  # 2
+        (0.0, 1.0, 0.0),  # 3
+    ]
+
+    el = [(0, 1, 2), (0, 2, 3)]
+
+    g = GeometryFM2D(node_coordinates=nc, element_table=el, projection="LONG/LAT")
+    g2 = GeometryFM2D(node_coordinates=nc, element_table=el, projection="LONG/LAT")
+
+    assert g == g2
+
+    g3 = GeometryFM2D(node_coordinates=nc, element_table=el, projection="UTM-33")
+    assert g != g3
+
+
+def test_equality_shifted_coords():
+    nc1 = [
+        (0.0, 0.0, 0.0),  # 0
+        (1.0, 0.0, 0.0),  # 1
+        (1.0, 1.0, 0.0),  # 2
+        (0.0, 1.0, 0.0),  # 3
+    ]
+
+    el = [(0, 1, 2), (0, 2, 3)]
+    g = GeometryFM2D(node_coordinates=nc1, element_table=el, projection="LONG/LAT")
+
+    nc2 = [
+        (0.1, 0.0, 0.0),  # 0
+        (1.0, 0.0, 0.0),  # 1
+        (1.0, 1.0, 0.0),  # 2
+        (0.1, 1.0, 0.0),  # 3
+    ]
+
+    g2 = GeometryFM2D(node_coordinates=nc2, element_table=el, projection="LONG/LAT")
+    assert g != g2