diff --git a/hydrolib/core/dflowfm/net/models.py b/hydrolib/core/dflowfm/net/models.py index 2e5836dcc..aea2a53fe 100644 --- a/hydrolib/core/dflowfm/net/models.py +++ b/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 diff --git a/hydrolib/core/dflowfm/net/reader.py b/hydrolib/core/dflowfm/net/reader.py index 45403b815..3a716922c 100644 --- a/hydrolib/core/dflowfm/net/reader.py +++ b/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: diff --git a/tests/dflowfm/test_net.py b/tests/dflowfm/test_net.py index 0792740d5..c2e175e4a 100644 --- a/tests/dflowfm/test_net.py +++ b/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():