Prevent duplicate contacts/links in Link1d2d

hydrolib/core/dflowfm/net/models.py

@@ -590,12 +590,29 @@ class Link1d2d(BaseModel):
 
  meshkernel: mk.MeshKernel = Field(default_factory=mk.MeshKernel)
 
- link1d2d_id: np.ndarray = Field(default_factory=lambda: np.empty(0, object))
- link1d2d_long_name: np.ndarray = Field(default_factory=lambda: np.empty(0, object))
- link1d2d_contact_type: np.ndarray = Field(
- default_factory=lambda: np.empty(0, np.int32)
- )
- link1d2d: np.ndarray = Field(default_factory=lambda: np.empty((0, 2), np.int32))
+ @property
+ def link1d2d(self) -> np.ndarray[np.int32]:
+ contacts = self.meshkernel.contacts_get()
+ link1d2d_arr = np.stack(
+ [contacts.mesh1d_indices, contacts.mesh2d_indices], axis=1
+ )
+ return link1d2d_arr
+
+ @property
+ def link1d2d_contact_type(self) -> np.ndarray[np.int32]:
+ contacts = self.meshkernel.contacts_get()
+ link1d2d_contact_type_arr = np.full(contacts.mesh1d_indices.size, 3)
+ return link1d2d_contact_type_arr
+
+ @property
+ def link1d2d_id(self) -> np.ndarray[object]:
+ link1d2d_id_arr = np.array([f"{n1d:d}_{f2d:d}" for n1d, f2d in self.link1d2d])
+ return link1d2d_id_arr
+
+ @property
+ def link1d2d_long_name(self) -> np.ndarray[object]:
+ link1d2d_id_arr = np.array([f"{n1d:d}_{f2d:d}" for n1d, f2d in self.link1d2d])
+ return link1d2d_id_arr
 
  def is_empty(self) -> bool:
  """Whether this Link1d2d is currently empty.
@@ -626,29 +643,6 @@ def clear(self) -> None:
  self.meshkernel._allocate_state(self.meshkernel.get_projection())
  self.meshkernel.contacts_get()
 
- def _process(self) -> None:
- """
- Get links from meshkernel and add to the array with link administration
- """
- contacts = self.meshkernel.contacts_get()
-
- self.link1d2d = np.append(
- self.link1d2d,
- np.stack([contacts.mesh1d_indices, contacts.mesh2d_indices], axis=1),
- axis=0,
- )
- self.link1d2d_contact_type = np.append(
- self.link1d2d_contact_type, np.full(contacts.mesh1d_indices.size, 3)
- )
- self.link1d2d_id = np.append(
- self.link1d2d_id,
- np.array([f"{n1d:d}_{f2d:d}" for n1d, f2d in self.link1d2d]),
- )
- self.link1d2d_long_name = np.append(
- self.link1d2d_long_name,
- np.array([f"{n1d:d}_{f2d:d}" for n1d, f2d in self.link1d2d]),
- )
-
  def _link_from_1d_to_2d(
  self, node_mask: np.ndarray, polygon: mk.GeometryList = None
  ):
@@ -668,7 +662,6 @@ def _link_from_1d_to_2d(
  self.meshkernel.contacts_compute_single(
  node_mask=node_mask, polygons=polygon, projection_factor=1.0
  )
- self._process()
 
  # Note that the function "contacts_compute_multiple" also computes the connections, but does not take into account
  # a bounding polygon or the end points of the 1d mesh.
@@ -678,7 +671,6 @@ def _link_from_2d_to_1d_embedded(
  ):
  """"""
  self.meshkernel.contacts_compute_with_points(node_mask=node_mask, points=points)
- self._process()
 
  def _link_from_2d_to_1d_lateral(
  self,
@@ -693,7 +685,6 @@ def _link_from_2d_to_1d_lateral(
  self.meshkernel.contacts_compute_boundary(
  node_mask=node_mask, polygons=polygon, search_radius=search_radius
  )
- self._process()
 
 
 class Mesh1d(BaseModel):

hydrolib/core/dflowfm/net/reader.py

@@ -121,12 +121,13 @@ def read_link1d2d(self, link1d2d: Link1d2d) -> None:
  ds = nc.Dataset(self._ncfile_path) # type: ignore[import]
 
  # Read mesh1d
- for meshkey, nckey in self._explorer.link1d2d_var_name_mapping.items():
-  setattr(link1d2d, meshkey, self._read_nc_attribute(ds[nckey]))
+ link1d2d_contact_type_arr = self._read_nc_attribute(ds["link1d2d_contact_type"])
+ assert (link1d2d_contact_type_arr == 3).all()
 
+ link1d2d_arr = self._read_nc_attribute(ds["link1d2d"])
  # set contacts on meshkernel, use .copy() to avoid strided arrays
- mesh1d_indices = link1d2d.link1d2d[:, 0].copy()
- mesh2d_indices = link1d2d.link1d2d[:, 1].copy()
+ mesh1d_indices = link1d2d_arr[:, 0].copy()
+ mesh2d_indices = link1d2d_arr[:, 1].copy()
  contacts = Contacts(
  mesh1d_indices=mesh1d_indices, mesh2d_indices=mesh2d_indices
  )

hydrolib/core/dflowfm/net/writer.py

@@ -330,6 +330,42 @@ def _set_1dmesh(self, ncfile: nc.Dataset, mesh1d: Mesh1d) -> None: # type: igno
  mesh1d_node_offset.units = "m"
  mesh1d_node_offset[:] = mesh1d.mesh1d_node_branch_offset
 
+ mesh_edge_x = ncfile.createVariable(
+ "mesh1d_edge_x", np.float64, "mesh1d_nEdges"
+ )
+ mesh_edge_x.standard_name = "projection_x_coordinate"
+ mesh_edge_x.long_name = (
+ "Characteristic x-coordinate of the mesh edge (e.g. midpoint)"
+ )
+ mesh_edge_x.units = "m"
+ mesh_edge_x[:] = mesh1d.mesh1d_edge_x
+
+ mesh_edge_y = ncfile.createVariable(
+ "mesh1d_edge_y", np.float64, "mesh1d_nEdges"
+ )
+ mesh_edge_y.standard_name = "projection_y_coordinate"
+ mesh_edge_y.long_name = (
+ "Characteristic y-coordinate of the mesh edge (e.g. midpoint)"
+ )
+ mesh_edge_y.units = "m"
+ mesh_edge_y[:] = mesh1d.mesh1d_edge_y
+
+ mesh_node_x = ncfile.createVariable(
+ "mesh1d_node_x", np.float64, "mesh1d_nNodes"
+ )
+ mesh_node_x.standard_name = "projection_x_coordinate"
+ mesh_node_x.long_name = "x coordinates of mesh nodes"
+ mesh_node_x.units = "m"
+ mesh_node_x[:] = mesh1d.mesh1d_node_x
+
+ mesh_node_y = ncfile.createVariable(
+ "mesh1d_node_y", np.float64, "mesh1d_nNodes"
+ )
+ mesh_node_y.standard_name = "projection_y_coordinate"
+ mesh_node_y.long_name = "y coordinates of mesh nodes"
+ mesh_node_y.units = "m"
+ mesh_node_y[:] = mesh1d.mesh1d_node_y
+
  def _set_2dmesh(self, ncfile: nc.Dataset, mesh2d: Mesh2d) -> None: # type: ignore[import]
 
  nc_mesh2d = ncfile.createVariable("mesh2d", "i4", ())

tests/dflowfm/test_net.py

@@ -98,9 +98,7 @@ def get_circle_gl(r, detail=100):
  return polygon
 
 
-@pytest.mark.plots
-def test_create_1d_2d_1d2d():
-
+def network_1d_2d_1d2dlinks():
  # Define line (spiral)
  theta = np.arange(0.1, 20, 0.01)
 
@@ -132,6 +130,12 @@ def test_create_1d_2d_1d2d():
 
  # Add links
  network.link1d2d_from_1d_to_2d(branchids=["branch1"], polygon=get_circle_gl(19))
+ return network
+
+
+@pytest.mark.plots
+def test_create_1d_2d_1d2d():
+ network = network_1d_2d_1d2dlinks()
 
  mesh2d_output = network._mesh2d.get_mesh2d()
  assert len(mesh2d_output.face_x) == 152
@@ -171,8 +175,25 @@ def test_create_1d_2d_1d2d():
  network2.plot(ax=ax2)
 
 
-def test_create_2d():
+def test_create_1d_2d_1d2d_call_link_generation_twice():
+ """
+ a second generation of links should overwrite the existing links,
+ this testcase checks whether that is the case.
+ Related issue: https://github.com/Deltares/HYDROLIB-core/issues/546
+ """
+ network = network_1d_2d_1d2dlinks()
 
+ # Add links again, do this twice to check if contacts are overwritten and not appended
+ network.link1d2d_from_1d_to_2d(branchids=["branch1"], polygon=get_circle_gl(19))
+
+ network1_link1d2d = network._link1d2d.link1d2d
+ assert network1_link1d2d.shape == (21, 2)
+ assert len(network._link1d2d.link1d2d_contact_type) == 21
+ assert len(network._link1d2d.link1d2d_id) == 21
+ assert len(network._link1d2d.link1d2d_long_name) == 21
+
+
+def test_create_2d():
  # Define polygon
  bbox = (1.0, -2.0, 3.0, 4.0)