reintroduced mesh2d_node_x and others as properties, fixes two testcases

hydrolib/core/dflowfm/net/models.py

@@ -121,7 +121,48 @@ class Mesh2d(
     # mesh2d_face_nodes: np.ndarray = Field(
     #     default_factory=lambda: np.empty((0, 0), dtype=np.int32)
     # )
-
+
+    @property
+    def mesh2d_node_x(self):
+        return self.meshkernel.mesh2d_get().node_x
+
+    @property
+    def mesh2d_node_y(self):
+        return self.meshkernel.mesh2d_get().node_y
+
+    @property
+    def mesh2d_edge_x(self):
+        return self.meshkernel.mesh2d_get().edge_x
+
+    @property
+    def mesh2d_edge_y(self):
+        return self.meshkernel.mesh2d_get().edge_y
+
+    @property
+    def mesh2d_edge_nodes(self):
+        mesh2d_output = self.meshkernel.mesh2d_get()
+        edge_nodes = mesh2d_output.edge_nodes.reshape((-1, 2))
+        return edge_nodes
+
+    @property
+    def mesh2d_face_x(self):
+        return self.meshkernel.mesh2d_get().face_x
+
+    @property
+    def mesh2d_face_y(self):
+        return self.meshkernel.mesh2d_get().face_y
+
+    @property
+    def mesh2d_face_nodes(self):
+        mesh2d_output = self.meshkernel.mesh2d_get()
+        npf = mesh2d_output.nodes_per_face
+        if self.is_empty():
+            return np.empty((0, 0), dtype=np.int32)
+        face_node_connectivity = np.full((len(mesh2d_output.face_x), max(npf)), np.iinfo(np.int32).min)
+        idx = (np.ones_like(face_node_connectivity) * np.arange(max(npf))[None, :]) < npf[:, None]
+        face_node_connectivity[idx] = mesh2d_output.face_nodes
+        return face_node_connectivity
+
     def is_empty(self) -> bool:
         """Determine whether this Mesh2d is empty.
 
@@ -140,23 +181,13 @@ def read_file(self, file_path: Path) -> None:
         reader = UgridReader(file_path)
         reader.read_mesh2d(self)
 
-    def _set_mesh2d(self, node_x=None, node_y=None, edge_nodes=None) -> None:
+    def _set_mesh2d(self, node_x, node_y, edge_nodes) -> None:
         # TODO: setting types is necessary since meshkernel.mesh2d_set requires them very specifically
-        #TODO: this func is only used once, while providing input, so consider making it obligatory arguments
-
-        if node_x is None:
-            node_x = self.mesh2d_node_x
-        if node_y is None:
-            node_y = self.mesh2d_node_y
-        if edge_nodes is None:
-            edge_nodes = self.mesh2d_edge_nodes.ravel()
-
         mesh2d = mk.Mesh2d(
             node_x=node_x.astype(np.float64),
             node_y=node_y.astype(np.float64),
             edge_nodes=edge_nodes.ravel().astype(np.int32),
         )
-
         self.meshkernel.mesh2d_set(mesh2d)
 
     def get_mesh2d(self) -> mk.Mesh2d:
@@ -236,21 +267,6 @@ def _process(self, mesh2d_input) -> None:  # TODO: input arg is not used, so rem
         # TODO: this is the mesh2d_face_node_connectivity, not mesh2d_face_nodes
         return
 
-    @property
-    def mesh2d_edge_nodes(self):
-        mesh2d_output = self.meshkernel.mesh2d_get()
-        edge_nodes = mesh2d_output.edge_nodes.reshape((-1, 2))
-        return edge_nodes
-
-    @property
-    def mesh2d_face_nodes(self):
-        mesh2d_output = self.meshkernel.mesh2d_get()
-        npf = mesh2d_output.nodes_per_face
-        face_node_connectivity = np.full((len(mesh2d_output.face_x), max(npf)), np.iinfo(np.int32).min)
-        idx = (np.ones_like(face_node_connectivity) * np.arange(max(npf))[None, :]) < npf[:, None]
-        face_node_connectivity[idx] = mesh2d_output.face_nodes
-        return face_node_connectivity
-
     def clip(
         self,
         geometrylist: mk.GeometryList,
@@ -763,8 +779,8 @@ class Mesh1d(BaseModel):
     network1d_edge_nodes: np.ndarray = Field(
         default_factory=lambda: np.empty((0, 2), np.int32)
     )
-    network1d_geom_x: np.ndarray = Field(default_factory=lambda: np.empty(0, np.double))
-    network1d_geom_y: np.ndarray = Field(default_factory=lambda: np.empty(0, np.double))
+    network1d_geom_x: np.ndarray = Field(default_factory=lambda: np.empty(0, np.double)) #TODO: sync with meshkernel.node_x
+    network1d_geom_y: np.ndarray = Field(default_factory=lambda: np.empty(0, np.double)) #TODO: sync with meshkernel.node_y
     network1d_part_node_count: np.ndarray = Field(
         default_factory=lambda: np.empty(0, np.int32)
     )
@@ -782,7 +798,7 @@ class Mesh1d(BaseModel):
         default_factory=lambda: np.empty(0, np.double)
     )
 
-    mesh1d_edge_nodes: np.ndarray = Field(
+    mesh1d_edge_nodes: np.ndarray = Field( #TODO: sync with meshkernel.edge_nodes (ravelled)
         default_factory=lambda: np.empty((0, 2), np.int32)
     )
     mesh1d_edge_x: np.ndarray = Field(default_factory=lambda: np.empty(0, np.double))
@@ -1255,6 +1271,16 @@ def mesh1d_add_branch(
         )
         self._mesh1d._set_mesh1d()
         return name
+
+    def plot(self):
+        import matplotlib.pyplot as plt
+        fig,ax = plt.subplots()
+        mesh2d_output = self._mesh2d.get_mesh2d()
+        mesh1d_output = self._mesh1d._get_mesh1d()
+        links_output = self._link1d2d.meshkernel.contacts_get()
+        mesh2d_output.plot_edges(ax=ax, color='r')
+        mesh1d_output.plot_edges(ax=ax, color='g')
+        links_output.plot_edges(ax=ax, mesh1d=mesh1d_output, mesh2d=mesh2d_output, color='k')
 
 
 class NetworkModel(ParsableFileModel):
@@ -1323,3 +1349,7 @@ def _get_serializer(cls):
     def _get_parser(cls):
         # Unused, but requires abstract implementation
         pass
+
+    @property
+    def plot(self):
+        return self.network.plot

hydrolib/core/dflowfm/net/reader.py

@@ -118,6 +118,15 @@ def read_link1d2d(self, link1d2d: Link1d2d) -> None:
         for meshkey, nckey in self._explorer.link1d2d_var_name_mapping.items():
             setattr(link1d2d, meshkey, self._read_nc_attribute(ds[nckey]))
 
+        # TODO: setting contacts is not possible in meshkernel (e.g. with contacts_set())
+        # so misalignment between link1d2d.link1d2d and 
+        # empty _link1d2d.meshkernel.contacts_get().mesh2d_indices
+        # mesh1d_indices = link1d2d.link1d2d[:,0]
+        # mesh2d_indices = link1d2d.link1d2d[:,1]
+        # import meshkernel as mk
+        # contacts = mk.Contacts(mesh1d_indices=mesh1d_indices, mesh2d_indices=mesh2d_indices)
+        # link1d2d.meshkernel.contacts_set(contacts)
+
         ds.close()
 
     def _read_nc_attribute(self, attr: nc._netCDF4.Variable) -> np.ndarray:

tests/dflowfm/test_net.py

@@ -17,15 +17,15 @@
 from ..utils import test_input_dir, test_output_dir
 
 
-def plot_network(network):
+def plot_network(network): #TODO: can be removed
     _, ax = plt.subplots()
     ax.set_aspect(1.0)
-    network.plot(ax=ax)  # TODO: does not exist
+    network.plot(ax=ax) #TODO: newly added but was already called
     ax.autoscale()
     plt.show()
 
 
-def _plot_mesh2d(mesh2d, ax=None, **kwargs):
+def _plot_mesh2d(mesh2d, ax=None, **kwargs): #TODO, this can be removed meshkernel.mesh2d_get().plot_edges() does the trick
     from matplotlib.collections import LineCollection
 
     if ax is None:
@@ -161,34 +161,36 @@ def test_create_1d_2d_1d2d():
     mesh1d_output = network._mesh1d._get_mesh1d()
     assert len(mesh1d_output.node_x) == 110
 
-    network_link1d2d = network._link1d2d.link1d2d
-    network_con_m1d = network._link1d2d.meshkernel.contacts_get().mesh1d_indices
-    network_con_m2d = network._link1d2d.meshkernel.contacts_get().mesh2d_indices
-    assert network_link1d2d.shape == (21, 2)
-    assert network_con_m1d.size == 21
-    assert network_con_m2d.size == 21
-
+    network1_link1d2d = network._link1d2d.link1d2d
+    network1_con_m1d = network._link1d2d.meshkernel.contacts_get().mesh1d_indices
+    network1_con_m2d = network._link1d2d.meshkernel.contacts_get().mesh2d_indices
+    assert network1_link1d2d.shape == (21, 2)
+    assert network1_con_m1d.size == 21
+    assert network1_con_m2d.size == 21
+    
     # Write to file
     file_out = Path(test_output_dir / "test_net.nc")
     network.to_file(file_out)
 
     #read from written file
     network2 = NetworkModel(file_out)
-
+    
     mesh2d_output = network2._mesh2d.get_mesh2d()
-    assert len(mesh2d_output.face_x) == 152
+    assert len(network2._mesh2d.mesh2d_face_x) == 152
     mesh1d_output = network2._mesh1d._get_mesh1d()
     assert len(mesh1d_output.node_x) == 110
 
-    network_link1d2d = network2._link1d2d.link1d2d
-    network_con_m1d = network2._link1d2d.meshkernel.contacts_get().mesh1d_indices
-    network_con_m2d = network2._link1d2d.meshkernel.contacts_get().mesh2d_indices
-    assert network_link1d2d.shape == (21, 2)
-    #TODO: below fail, so the meshkernel instance is not up to date for some reason, should it be?
-    assert network_con_m1d.size == 21
-    assert network_con_m2d.size == 21
+    network2_link1d2d = network2._link1d2d.link1d2d
+    network2_con_m1d = network2._link1d2d.meshkernel.contacts_get().mesh1d_indices
+    network2_con_m2d = network2._link1d2d.meshkernel.contacts_get().mesh2d_indices
+    assert network2_link1d2d.shape == (21, 2)
+    #TODO: below asserts fail, since the meshkernel contacts are not set upon reading (not implemented in meshkernel)
+    assert network2_con_m1d.size == 21
+    assert network2_con_m2d.size == 21
 
-    # plot_network(network)
+    # plot both networks
+    network.plot()
+    network2.plot()
 
 
 def test_create_2d():